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14.  Abstract 

This  Workshop  focused  on  strategies  for  promoting  and  developing  engineering-level  transonic 
flutter  prediction  techniques. 

The  technology  of  transonic  aerodynamics  is  currently  undergoing  rapid  development.  Significant 
progress  is  being  made  to  solve  the  inherently  nonlinear  equations  describing  unsteady  motions  of 
wings  in  transonic  flow,  while  the  availability  of  reliable  and  efficient  computational  methods  will 
greatly  enhance  the  ability  to  predict  the  aeroelastic  behaviour  of  modem  aircraft  operating  under 
transonic  flow  conditions. 

AGARD-SMP  has  previously  coordinated  unsteady  aerodynamic  research  carried  out  on  a 
number  of  “standard”  wind  tunnel  model  configurations  and  published  the  results.  The  proposals 
contained  in  the  Evaluation  Report  (WJ.Mykytow)  on  the  Fall  1984  Structures  and  Materials 
Panel  Specialists’  Meeting  on  “Transonic  Unsteady  Aerodynamics”,  together  with  an  expanded 
range  of  aeroelastic  configurations,  formed  the  guidelines  which  this  Workshop  followed. 

Two  of  the  keynote  papers  presented  at  this  Workshop  are  published  in  this  Report. 
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PREFACE 


The  technology  of  transonic  aerodynamics  is  currently  undergoing  rapid  development.  Significant  progress  is  being  made 
to  solve  the  inherently  nonlinear  equations  describing  unsteady  motions  of  wings  in  transonic  flow,  while  the  availability  of 
reliable  and  efficient  computational  methods  will  greatly  enhance  the  ability  to  predict  the  aeroelastic  behaviour  of  modem 
aircraft  operating  under  transonic  flow  conditions. 

Calculation  of  the  motion-induced  unsteady  airloads  on  oscillating  wings  in  transonic  flow  requires  large  computer 
storage  capacity  and  long  computing  times.  Several  computation  techniques  involving  various  degrees  of  approximation  have 
been  elaborated  and  are  now  available  for  use.  However,  it  is  difficult  to  compare  and  evaluate  all  these  competing  methods 
with  regard  to  their  applicability  for  economic,  engineering-type  aeroelastic  predictions.  Indeed,  from  the  relatively  simple 
Transonic  Small  Perturbation  (TSP)  equation  up  to  the  Full  Potential  Euler  and  Navier-Stokes  Equations,  a  wide  range  of 
unsteady  transonic  flow  prediction  capabilities  exists  which  apply  finite  difference,  finite  element  and  integral  equation 
techniques,  the  results  of  which  are  manifested  not  only  in  computational  accuracy  but  also  in  large  differences  of  computing 
time  and  computer  requirements.  This  was  clearly  demonstrated  at  the  first  AGARD-SMP  Specialists’  Meeting  on  “Unsteady 
Airloads  in  Separated  and  Transonic  Flow”  held  in  Lisbon,  Portugal,  in  April  1977.  At  its  subsequent  fall  1977  meeting,  the 
AGARD-SMP  then  formed  a  working  group  on  “Standard  Configurations  for  Aeroelastic  Applications  of  Transonic 
Unsteady  Aerodynamics”  with  members  of  six  NATO  countries  under  the  auspices  of  the  Subcommittee  on  Aero  elasticity. 
The  aim  of  this  working  group  was  to  accelerate  —  in  a  cooperative  program  —  the  development  of  new  theoretical,  numerical 
and  experimental  techniques  in  transonic  unsteady  aerodynamics  and  their  application  to  aeroelastic  problems  of  aircraft 
loads,  stability  and  flutter.  Based  on  the  recommendations  obtained  by  the  members  of  this  working  group  from  aeroelasticians 
and  aerodynamicists  in  their  respective  countries,  a  limited  set  of  two-dimensional  and  three-dimensional  aeroelastic 
configurations  was  defined  as  standard  airfoils  and  wings,  and  a  compendium  of  unsteady  aerodynamic  measurements  in  the 
transonic  flow  regime  was  elaborated. 

After  AGARD  had  established  these  standard  configurations,  the  next  step  was  to  request  aeroelasticians  in  the  NATO 
countries  to  evaluate  their  methods.  This  effort  culminated  in  a  Specialists’  Meeting  on  ‘Transonic  Unsteady  Aerodynamics 
and  Its  Aeroelastic  Applications”  held  at  the  SMP  fall  1984  meeting  in  Toulouse.  This  meeting  clearly  demonstrated  that  the 
AGARD  cooperative  program  on  transonic  unsteady  aerodynamics  and  aeroelastic  applications  has  indeed  stimulated 
excellent  contributions  to  advance  the  state  of  the  art  which  has  progressed  substantially  in  the  seven  years  between  the  last  two 
SMP  Specialists’  Meetings  in  spring  1977  and  fall  1984.  Moreover,  it  was  also  shown  at  the  last  meeting  that  encouraging  three- 
dimensional  methods  are  now  emerging  to  predict  unsteady  airloads  for  transonic  flutter  analyses  of  clean  wings.  In  order  to 
assess  the  applicability  of  these  prediction  techniques  for  engineering-type  practical  aeroelastic  analysis,  the  desirability  of  re¬ 
establishing  a  set  of  standard  configurations  for  comparisons  of  calculated  and  measured  dynamic  aeroelastic  behaviour  under 
high  subsonic  to  transonic  flow  conditions  was  discussed  at  the  fall  1984  meeting  of  the  AGARD-SMP  Subcommittee  on 
Aeroelasticity.  Hence,  two  2-D  airfoils  were  chosen  in  an  initial  selection  at  the  last  SMP  meeting  in  spring  1986  as  transonic 
aeroelastic  standard  configurations.  Assuming  further  progress  in  the  development  of  transonic  flutter  prediction  methods  it  is 
hoped  that  the  selection  of  these  standard  configurations  will  be  completed  by  including  a  more  straightforward  3-D  wing  with 
supercritical  profile,  a  low  aspect  ratio  fighter- type  swept  wing,  a  wing  with  controls,  etc.  The  cooperative  AGARD  program 
and  activities  based  on  and  dictated  by  these  standard  aeroelastic  configurations,  then,  should  culminate  in  another  symposium 
on  ‘Transonic  Unsteady  Aerodynamics  and  Its  Aeroelastic  Applications”  to  be  held  in  1990—1991. 

In  order  to  discuss  strategies  for  promoting  engineering-level  transonic  aeroelastic  prediction  techniques  and  to  establish 
tops  for  further  research  work,  a  Workshop  entitled 

“Future  Research  on  Transonic  Unsteady  Aerodynamics 
and  Its  Aeroelastic  Applications” 

has  been  held  under  the  auspices  of  the  SMP- Subcommittee  on  Aeroelasticity  at  the  fall  1986  AGARD-SMP  Meeting  in 
Athens.  Two  of  the  keynote  papers  presented  at  this  Workshop  are  published  in  this  Report. 


H.Forsching 

Chairman,  Subcommittee  on 
Aeroelasticity 


ABSTRACT 


This  Workshop  focused  on  strategies  for  promoting  and  developing  engineering-level  transonic  flutter  prediction 
techniques. 

The  technology  of  transonic  aerodynamics  is  currently  undergoing  rapid  development.  Significant  progress  is  being  made 
to  solve  the  inherently  nonlinear  equations  describing  unsteady  motions  of  wings  in  transonic  flow,  while  the  availability  of 
reliable  and  efficient  computational  methods  will  greatly  enhance  the  ability  to  predict  the  aeroelastic  behaviour  of  modem 
aircraft  operating  under  transonic  flow  conditions. 

AGARD-SMP  has  previously  coordinated  unsteady  aerodynamic  research  carried  out  on  a  number  of  “standard”  wind 
tunnel  model  configurations  and  published  the  results.  The  proposals  contained  in  the  Evaluation  Report  (W J.Mykytow)  on 
the  Fall  1984  Structures  and  Materials  Panel  Specialists’  Meeting  on  ‘Transonic  Unsteady  Aerodynamics”,  together  with  an 
expanded  range  of  aeroelastic  configurations,  formed  the  guidelines  which  this  Workshop  followed. 

Two  of  the  keynote  papers  presented  at  this  Workshop  are  published  in  this  Report. 


*  *  * 


Cet  Atelier  a  concentre  ses  efforts  sur  les  strategies  visant  a  promouvoir  et  a  developper  des  techniques  de  prevision  du 
flottement  en  vol  transsonique  au  niveau  de  Fingenierie. 

La  technologie  de  Faerodynamique  transsonique  fait  actuellement  l’objet  d’un  developpement  rapide.  Des  progres 
significatifs  sont  realises  dans  la  resolution  des  equations  fondamentalement  nonlineaires  decrivant  les  mouvements 
instationnaires  des  ailes  dans  le  flux  transsonique,  alors  que  la  disponibilite  de  methodes  de  calcul  fiables  et  efficaces 
augmenteront  considerablement  la  possibility  de  prevoir  le  comportement  aeroelastique  des  avions  modemes  volant  dans  des 
conditions  de  flux  transsonique. 

Le  Groupe  charge  des  Structures  et  des  Materiaux  de  l’AGARD  a  precedemment  coordonne  les  travaux  de  recherche  sur 
Faerodynamique  instationnaire  effectues  sur  un  certain  nombre  de  configurations  de  maquettes  de  soufflerie  “standards”  et  en 
a  publie  les  resultats.  Les  propositions  contenues  dans  le  rapport  devaluation  (W  J.Mykytow)  etabli  a  Tissue  de  la  reunion,  a 
Fautomne  1984,  des  specialistes  du  Groupe  charge  de  l’etude  des  Structures  et  des  Materiaux,  sur  TAerodynamique 
Transsonique  Instationnaire”,  ainsi  qu’une  vaste  gamme  de  configurations  aeroelastiques,  ont  ete  successivement  examinees 
par  FAtelier  au  cours  de  ses  travaux. 

Deux  des  documents  clefs  presentes  au  cours  de  la  reunion  de  cet  Atelier  sont  publies  dans  ce  Rapport. 
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SUMMARY 

The  current  scope,  recent  progress,  and  plans  for  research  and  development  of 
computational  methods  for  unsteady  aerodynamics  at  the  NASA  Langley  Research  Center  are 
reviewed.  Both  i ntegral -equati on  and  finite-difference  methods  for  inviscid  and 
viscous  flows  are  discussed.  Although  the  great  bulk  of  the  effort  has  focused  on 
finite-difference  solution  of  the  transonic  smal 1  - perturbati on  equation,  the 
integral-equation  program  is  given  primary  emphasis  here  because  it  is  less  well  known. 


INTRODUCTION 

Progress  in  the  development  of  computational  methods  for  steady  and  unsteady 
aerodynamics  has  perennially  paced  advancements  in  aeroelastic  analysis  and  design 
capabilities.  These  capabilities,  in  turn,  are  of  growing  importance  in  the  analysis 
and  design  of  high-performance  aircraft  as  well  as  other  types  of  flight  vehicles. 
Consequently,  considerable  effort  has  been  directed  toward  the  development  of 
appropriate  unsteady- aerodynami c  methodology  in  the  NATO  countries  and  elsewhere.  This 
paper  reviews  the  contributions  to  those  efforts  at  the  NASA  Langley  Research  Center. 
Specifically,  the  current  scope,  recent  progress,  and  plans  for  research  and 
development  of  both  i ntegral - equati on  and  f i ni te- di f f erence  methods  for  inviscid  and 
viscous  flows  are  discussed,  and  example  applications  are  shown.  Although  the  great 
bulk  of  the  effort  in  recent  years  has  focused  on  f i ni te- di ff erence  solution  of  the 
transonic  smal 1  - perturbati on  equation,  the  i ntegral - equati on  program  is  given  primary 
emphasis  here  because  it  is  less  well  known. 


INTEGRAL-EQUATION  METHODS 

The  Langley  i ntegral - equati ons  program  is  directed  toward  general,  accurate,  efficient, 
and  unified  treatment  of  flows  around  vehicles  having  arbitrary  shapes,  motions,  and 
deformations  (including  control  motions)  at  subsonic,  transonic,  and  supersonic  speeds 
up  to  high  angles  of  attack.  Special  attention  is  given  to  real-world  design  and 
operating  conditions  (e.g.,  Mach  number,  angle  of  attack,  maneuver)  as  well  as  to 
efficient  computation  for  both  design  and  analysis  applications.  As  will  be  brought 
out  in  the  subsequent  discussion,  the  i ntegral - equati on  approach  is  well  suited  for 
these  purposes  because  flow  complexities  such  as  viscous  effects  or  transonic  flow  need 
to  be  addressed  only  in  the  flow  regions  where  they  actually  occur,  and  there  is  no 
requirement  for  patching  and  matching  flow  domains  or  regional  solutions.  Moreover, 
for  design  applications  repetitive  and  non repeti ti ve  portions  of  the  computations  are 
readily  separable,  and  the  required  sensitivities  of  aerodynamic  parameters  to 
variations  in  aircraft  geometry  can  be  readily  calculated.  Although  the 
i ntegral -equati ons  research  program  has  been  given  only  limited  and  intermittent 
support  for  the  last  several  years,  it  has  nevertheless  produced  some  significant 
results. 

Following  a  long-range  plan  established  a  number  of  years  ago  (fig.  1),  initial  efforts 
addressed  the  development  of  surf ace- panel  methods  for  subsonic  (refs.  1  to  5)  and 
supersonic  (refs.  1,  2,  6,  7)  linearized  potential  flow.  Current  activities  include 
nonlinear  methods  implementing  the  f ul 1  - potenti al  equation  for  hi gh- subsoni c/ transoni c/ 
low-supersonic  speeds  (ref.  8).  Although  the  initial  hi gh- subsoni c/ transoni c 
proof-of-concept  codes  (refs.  9  to  11)  implemented  the  smal 1  - perturbati on  potential 
equation,  there  is  noparticular  benefit  in  refining  codes  for  small-perturbation 
conditions  or  two- dimi nsi onal  flow  as  stepping-stones  toward  more  realistic 
conditions.  Consequently,  these  items  have  been  deleted  from  the  original  plan  (fig. 
1).  Some  computations  for  two-dimensional  flow  are  made  in  order  to  conserve  computer 
resources,  however. 
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Another  change  from  the  original  plan  (fig.  1)  shows  that  the  Euler  equations  are  not 
addressed  explicitly  in  this  program.  Modification  of  the  f ul 1  - potenti al  equation  to 
account  for  entropy  changes  across  shock  waves  (e.g.,  as  in  ref.  12)  should  greatly 
expand  the  usefulness  of  potenti al - fl ow  solutions  well  into  the  range  of  flow 
conditions  that  would  otherwise  require  Euler  solutions.  Consequently,  for  present 
purposes  we  go  directly  from  the  modified  ful 1  - potenti al  method  to  the  Navi er-Stokes 
equations  which  are  being  addressed  by  use  of  the  classical  Helmholtz 
scalar/vector-potential  decomposition  (refs.  13  to  15).  Euler  solutions  may,  of 
course,  be  obtained  from  Navi er-Stokes  methods  with  zero  viscosity.  We  are 
specifically  concerned  with  several  types  of  viscous  influences:  Thin  wakes  separating 
from  1 i f ti ng- su rf ace  edges  are  well  represented  by  i nvi sci d- f 1 ow  singularities  (vortex 
sheets).  Other  viscous  influences  require  solution  of  N a vi er- S tok es  equations  or 
equivalent.  These  influences  include  boundary-1 ayer  effects,  especially  on  deflected 
and/or  deflecting  control  surfaces,  shock/boundary- 1  ay er  interaction,  and  large  areas 
of  flow  separation  in  general.  Specifics  of  these  problem  areas  are  addressed  in  the 
subsequent  sections  of  this  paper. 


The  time-dependent  full - potenti al  partial  differential  equation 
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is  the  governing  equation  for  most  of  the  work  described  herein.  Application  of  the 
generalized  Green ' s- functi on  method  to  this  equation  yields  an  equivalent  integral 
equation  for  the  velocity  potential  <f>  at  any  point  P  in  the  flow  or  on  the  surface  of  a 
body  in  the  flow  at  any  time  t  (ref.  1). 
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where  <f>  is  the  perturbation  velocity  potential,  G  is  the  Green's  function,  F  represents 
all  the  nonlinear  terms,  a*,  is  the  freestream  speed  of  sound,  U*,  is  the  freestream 
speed,  x  is  the  coordinate  in  the  freestream  direction,  S(x,y,z,t)  =  0  defines  the  body 
surface,  and 

I  □  s|  ■  7  sx  *  S5  *  sz  4  s? 

The  exact  boundary  condition  on  the  body  surface  is 


The  time  integration  with  respect  to  ti  in  equation  (2)  is  made  trivial  by  choice  of 
a  subsonic  or  supersonic  source  pulse  as  the  Green's  function. 

An  important  point  here  is  that  only  the  nonlinear  terms  need  to  be  integrated  over  a 
fluid  volume.  The  linear  terms  are  integrated  only  over  the  surface  of  the  body  and 
its  wake.  Note  also  that  the  Green's  function  is  a  function  of  freestream  Mach  number, 
not  local  Mach  number.  Equations  (2)  and  (3)  have  been  formulated  and  computationally 
implemented  in  a  moving  frame  of  reference  so  that  they  are  applicable  to  problems  such 
as  helicopter  rotors  and  maneuvering  aircraft  as  well  as  aircraft  in  uniform  motion 
( refs.  16  to  19) . 


Linearized  Theory 


If  perturbations  from  freestream  velocity  are  small,  and  Mach  number  is  not  near  one 
nor  too  high  in  the  supersonic  range,  the  non-linear  terms  are  negligible,  and  the 
volume  integral  can  be  ignored.  The  remaining  surface  integral  of  the  linear  terms  is 
discretized  by  surface  paneling  (ref.  2)  (e.g.,  arbitrary  twisted  quadri 1 ateral  panels 
as  in  refs.  3  and  4).  The  unsteady -flow  solution  can  then  be  obtained  directly  by 
integration  in  time  domain,  or  a  time  solution  by  Laplace  transform  (refs.  2  to  4) 

converts  to  a  compl ex- frequency  domain  formulation  which  is  generally  more  efficient 
for  use  in  solving  linear  aeroelastic  problems. 


The  velocity  potential  on  the  paneled  surface  is  then  found  in  terms  of  the  normalwash 
distribution  which,  in  general,  is  known  from  the  input  shape,  ori entati on ,moti on ,  and 
deformation. 
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6jh  is  Kronecker  delta,  s  is  the  Laplace  transform  variable  (complex  frequency). 


^jh»Cjh»Djh,Fin,Gjn  are  integrals  over  surface  panels,  ejh,njn  are  la9 
functions,  ana  S^^-  +_  1  for  panels  adjacent  to  a  trailing  edge  on  upper  or  lower 


surface  of  the  body  and  is  zero  otherwise, 
potential  by  use  of  Bernoulli's  equation. 


Surface  pressures  are  obtained  from  the 


Several  features  of  equations  (4)  to  (6)  are  significant.  First,  the  elements  of  the  Y 
and  2  influence  matrices  are  independent  of  the  normalwash  and  hence  independent  of  the 
mode  of  motion  or  deflection.  Moreover,  these  matrix  elements  are  simple  functions  of 
the  complex  frequency  s  so  that  the  cost  of  changing  frequency  or  calculating  for 

multiple  frequencies  is  small.  The  influence  integrals  B,  C,  and  D  represent  integrals 

of  source,  doublet,  and  "ratelet"  distributions  over  each  body-surface  panel,  and 
integrals  F  and  G  are  the  corresponding  doublet  and  "ratelet"  integrals  for  wake 
panels.  For  a  given  paneling  geometry,  all  of  these  integrals  are  functions  only  of 
Mach  number.  If  a  problem  (e.g.,  dynamic  response  or  flutter)  involves  multiple  modes 
of  normalwash,  the  normalwash  vector  in  the  equation  becomes  a  matrix  of  modal  columns, 
and  the  potential  distributions  for  all  the  modes  can  be  found  in  a  single  solution. 

Similarly,  solutions  for  additional  modes  or  revised  modes  (as  in  a  struc tu ral - desi gn 

optimization  problem)  can  be  obtained  without  recalculating  the  Y  and  Z  matrices.  For 
use  in  design  processes,  this  formulation  also  appears  to  provide  a  general  and  very 
efficient  means  for  evaluating  sensitivities,  i.e.,  changes  in  aerodynamic  properties 
caused  by  changes  in  external  shape.  Demonstration  calculations  have  been  initialed. 

The  generality  and  versatility  of  this  approach  is  indicated  by  its  use  by  Rockwell 
International  for  flutter  analysis  of  the  space  shuttle  (fig. 2)  in  the  mid  1970's. 
Nearly  800  panels  were  used  on  the  orbiter,  and  up  to  60  modes  of  motion  were  used  in 
both  symmetric  and  antisymmetric  flutter  analyses.  Subsequently,  the  external  tank  and 
solid  rocket  boosters  were  added,  and  the  calculations  were  repeated  for  the  entire 
launch  configuration. 


For  development  purposes  equations  (4)  to  (6)  have  been  implemented  in  a  prototype  code 
called  SOUSSA  Pl.l  (Steady,  Oscillatory,  and  Unsteady  Subsonic  and  Supersonic 
Aerodynamics  -  Version  1.1)  Trefs.  3  and  4)  wFich  is  applicable  to  vehicles  having 
arbitrary  shapes,  motions,  and  deformations  in  subsonic  flow  only.  The  Pl.l  code 
employs  zeroth-order .( constant-potential )  panels  along  with  the  data  base  and 
da ta- handl 1 ng  utilities  of  the  SPAR  f i ni te- el ement  structural - analy si s  program.  These 
were  incorporated  because  SOUSSA  Pl.l  originally  was  intended  for  the  calculation  of 
steady-state  structural  loads  and  unsteady  aerodynamics  for  flutter  and  gust- response 
calculation  in  multidisciplinary  structural -optimi zati on  computations  employing  the 
SPAR  structural  analysis.  The  SPAR  components,  however,  are  unnecessary  for 
stand-alone  use.  More  efficient  data  handling  methods  for  stand-alone  operation  are 
available. 


Subsequent  to  the  completion  of  SOUSSA  Pl.l  several  significant  improvements  have  been 
incorporated,  and  others  have  been  defined  (ref.  5).  Among  the  latter  are 
implementation  of  higher-order  panels,  elimination  of  the  SPAR  components, 
transposition  and  revision  of  the  solution  algorithm  to  substantially  reduce 
input/output  operations,  and  improved  implementation  of  the  trai 1 i ng-edge  (Kutta)  flow 
condition. 

Some  program  improvements  already  incorporated  in  the  SOUSSA  code  include  the 
development  of  an  "out-of-core"  solver  to  permit  the  use  of  paneling  schemes  that  lead 
to  coefficient  matrices  too  large  to  fit  in  the  memory  of  modest-size  computers;  the 
replacement  of  the  paneled  wake  by  an  analytical  wake  (reducing  the  cost  of  a  typical 
run  by  about  one-half)  but  retaining  an  option  to  use  paneled  wakes  if  needed  (e.g., 
when  there  is  another  lifting  surface  in  the  wake);  and  replacing  the  rectangular 
integration  of  pressures  by  a  Gaussian  guadrature  scheme  to  improve  the  accuracy  of  the 
calculated  generalized  aerodynamic  forces.  These  improvements  are  incorporated  in  a 
replacement  for  the  SOUSSA  code  (called  UTSA)  which  is  under  development  at  a  low  level 
of  effort. 


Figure  3  (reproduced  from  ref.  5)  compares  a  chordwise  distribution  of  pressure 
coefficient  Cp  calculated  by  the  SOUSSA  surface-panel  method  with  pressures  measured 
on  a  clipped  delta  wing  oscillating  in  pitch  (ref.  20).  The  wing  had  a 
six-percent- thick  circular-arc  airfoil.  The  agreement  is  good  and  is  representative  of 
results  obtained  with  this  code.  Figure  4  compares  calculated  and  measured  steady 
upper-surface  pressures  at  two  chordwise  locations  x  on  an  outboard  station  (y=0.85)  on 
the  same  clipped  delta  wing.  Two  points  are  to  be  made:  First,  in  the  range  of  angle 
of  attack  a  (-2  deg  to  +2  deg)  where  pressure  varies  linearly,  the  agreement  is 
excellent.  Second,  for  this  sharp-edge  wing,  the  influence  of  the  leading-edge  vortex 
is  substantial  and  begins  at  a  low  angle  of  attack.  The  latter  behavior  emphasizes  the 
importance  of  our  treatment  of  vortex-type  flow  separation  to  be  discussed  below.  A 
phenomenological  description  of  the  relation  between  the  vortex  development  and  the 
pressure  variation  shown  is  given  in  ref.  21  (from  which  figure  4  was  taken)  and  in 
Appendix  A  of  ref.  20. 

In  addition  to  the  subsonic  capability  of  the  SOUSSA  program,  a  supersonic 
proof-of-concept  surface-panel  code  has  been  written  to  implement  1 i near- theory 
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algorithms  developed  in  refs.  6  and  7.  The  code  employs  first-order  panels  and,  like 
SOUSSA,  is  applicable  to  vehicles  having  arbitrary  shapes,  motions,  and  deformations. 
Validation  and  application  of  the  code  have  begun. 

The  only  significant  difference  between  subsonic  and  supersonic  formulations  is  in  the 
expressions  for  the  influence  integrals  B,C,D,  F,G  in  equations  (5)  and  (6)  (see,  e.g., 
ref. 2).  Other  portions  of  the  computations,  such  as  paneling  geometry  and  solution 
alorithms  are  common  to  both.  Consequently,  it  is  possible  that  the  computational 
capability  for  supersonic  flow  derived  from  this  proof-of-concept  code  will 
subsequently  be  incorporated  into  the  subsonic  code  UTSA. 

The  status  and  near-term  plans  for  1 i near- theory  surf ace- panel  methods,  which  are 
applicable  to  vehicles  having  arbitrary  shapes,  motions,  and  deformations,  may  be 
summarized  as  follows:  As  planned  the  SOUSSA  program  will  be  superseded  by  an  improved 
program  UTSA  which  incorporates  first-order  panels  as  well  as  other  improvements 
indicated  by  earier  work  with  SOUSSA.  Ultimately,  the  code  may  include  both  subsonic 
and  supersonic  capabilities.  Frequency-domain  computations  are  most  efficient  for 
implementing  linear  theory,  but  a  time-domain  version  is  also  retained  for  evaluation 
of  the  surface  integral  in  the  nonlinear  methods  described  next.  Specific  activities 
include  configuring  the  UTSA  code  for  efficient  use  in  interdisciplinary  design 
processes,  incorporating  special  elements  to  improve  accuracy  and  efficiency  near 
normalwash  discontinuities  (e.g.,  at  control  surfaces),  completing  the  initial 
demonstration  of  the  efficient  computation  of  sensitivities  of  aerodynamic  pressures 
and  loads  to  variations  in  planform,  and  general  check  out  and  validation. 

Nonl i near  Theory 

When  the  flow  approaches  transonic  conditions  and/or  flow  perturbations  (e.g.,  angle  of 
attack)  become  large,  the  nonlinear  terms  represented  by  F  in  equation  (2)  are  no 
longer  negligible,  and  the  volume  integral  must  be  evaluated  in  combination  with  the 
su rf ace- panel  evaluation  of  the  linear  terms  (refs.  9  to  11).  For  nonlinear  problems 
it  is  important  to  note  (1)  that  the  Green’s  function  depends  on  freestream  Mach 
number,  not  local  Mach  number,  and  (2)  that  the  integrand  of  the  volume  integral 
diminishes  rapidly  in  magnitude  with  increasing  distance  from  the  body  and  its  wake. 

For  application  to  nonlinear  problems  the  integral-equation  method  has  several  features 
which  make  it  particularly  attractive  for  general,  efficient  computational 
implementation:  (1)  Evaluation  of  an  integral  is  required  rather  than  the  numerical 

solution  of  a  partial  differential  equation,  which  is  a  more  sensitive  process.  (2) 

The  volume  integral  need  be  treated  only  in  the  limited  region  of  flow  in  which 
nonlinear  terms  are  of  significant  magnitude  rather  than  over  an  entire  computational 
domain.  In  fact,  as  the  integration  proceeds  away  from  the  body,  it  is  terminated  when 
the  integrand  falls  below  a  preselected  threshhold  value.  (3)  Required  accuracy  can  be 
attained  with  relatively  few  computational  grid  points  in  the  fluid  (computational 
domain  of  the  volume  integral).  (4)  The  code  is  numerically  stable  even  when 
moderate- to- 1 arge  time  steps  are  employed.  (5)  Correct  far-field  boundary  conditions 
are  automatically  imposed.  This  condition  is  particularly  important  for  unsteady 
flow.  Li  near- theory  behavior  in  the  far  field  is  inherent  in  the  integral-equation 
solution.  (6)  When  viscous  flows  are  treated  by  the  seal ar/vector- potenti al 
decomposition  (to  be  discussed  below),  interfacing  (patching  and  matching)  of  regional 
solutions  (e.g.,  inner  viscous  solution  and  outer  inviscid  solution)  is  not  required. 

(7)  Even  for  solution  of  the  ful 1  - potenti al  equation,  there  is  no  requirement  for 
generating,  imbedding,  or  interpolating  surf ace- fi tted  computational  grids. 

In  this  section  smal 1 -perturbati on  transonic  attached  flow  will  be  considered  first 
followed  by  1 arge- perturbati on  subsonic  and  transonic  flow  conditions  involving 
vortex-type  flow  separation  in  the  form  of  thin  wakes  emanating  from  1 i f ti ng- su rf ace 
edges  and  finally  flow  conditions  involving  significant  viscous  effects  which  require 
solution  of  Navi er-Stokes  equations  for  attached  or  separated  flow  for  which  the 
seal ar/vector- potenti al  method  is  employed. 

Smal 1 -P erturbati on  Transonic  Flow:  For  proof-of-concept  demonstration  of  transonic 
capability,  only  the  TinaTT - perturbati  on  terms  were  retained  in  the  volume  integral  of 
equation  (2),  and  the  resulting  time-domain  computer  code  (ref.  11)  was  called  SUSAN 
(Steady  and  Unsteady  Subsonic  Aerodynamics-Nonlinear).  Figure  5  shows  chordwise 
pressure  distribution  near  the  root  of  a  rectangular  wing  as  calculated  by  the  SUSAN 
code  and  by  a  transonic  small-perturbation  f i ni te- di f f ere nee  code.  The  shock  is 
captured,  and  the  agreement  is  quite  good  even  though  only  a  few  elements  were  used  to 
evaluate  the  volume  integral,  and  the  domain  of  integration  extended  only  one  chord 
length  from  the  wing  perimeter.  Good  agreement  with  measured  pressures  (from  ref.  22) 
is  shown  in  figure  6  for  a  sharp-edge  wing  under  conditions  involving  supercritical 
flow  over  much  of  the  chord. 

Evolution  of  the  lifting  pressure  ACp  on  a  wing  oscillating  slowly  in  pitch  about  the 
leading  edge  is  shown  in  figure  7  at  three  times  during  a  cycle  of  motion.  Although 
only  ten  computational  elements  along  the  wing  chord  were  used  to  evaluate  the  volume 
integral  of  the  nonlinear  terms,  the  build-up  of  lift  and  the  appearance  of  a  shockwave 
are  clearly  indicated.  In  this  particular  figure,  the  symbols  shown  are  used  only  to 
distinguish  the  curves  and  do  not  indicate  computational  points. 
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The  formulation  described  here  and  its  implementation  in  the  SUSAN  code  demonstrated 
the  merits  of  the  i ntegral -equati on  method  for  transonic  flow.  However,  no  further 
development  of  the  small -perturbati on  approximation  is  planned. 


Subsonic/Transonic  Flow  with  Vortex  Separation:  All  of  the  preceding  involved 
calculation  of  the  velocity  potential.  For  solving  nonlinear  problems,  however,  there 
are  advantages  in  calculating  velocities  directly,  especially  when  large  velocity 
variations  occur,  when  shocks  are  present,  when  thin-wake  (vortex-like)  flow  separation 
from  wing  leading  or  side  edges  occurs  (fig.  8),  or  even  when  trai 1 i ng-edge  wake 
deformations  are  significant.  Taking  the  gradient  of  the  integral  equation  for  the 
potential  (equation  (2))  or  alternatively  applying  the  Green 1 s- functi on  method  to  the 
f ull -potenti al  equation  in  the  form  (for  steady  state) 


gives  (ref.  8) 
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where  p  is  the  fluid  density,  u>  is  the  vorticity  vector,  R  is  the  vector  from  "sending" 
point  to  "receiving"  point,  Er  is  a  unit  vector  in  the  R  direction,  and 
Eoo  is  a  unit  vector  in  freestream  direction. 

Equation  (8)  is  an  expression  for  the  velocity  field  V  as  the  sum  of  four  components: 
(1)  freestream,  (2)  a  surface  integral  which  gives  the  velocity  induced  by  the  flow 
singularities  representing  the  solid  body,  (3)  a  surface  integral  which  gives  the 
velocity  induced  by  the  vorticity  representing  the  thin  wake,  and  (4)  a  volume  integral 
representing  the  compressibility  terms  (right-hand  side  of  equation  (7)).  The 
integrand  of  this  volume  integral  decreases  more  rapidly  than  the  square  of  the 
distance  from  the  body  or  vortex  surface,  so  the  domain  of  integration  can  be 
relatively  small.  The  integrands  in  the  three  integrals  are  not  independent,  and 
solution  is  by  iteration  to  satisfy  the  boundary  conditions  on  the  body  and  to  deform 
the  free  vortex  sheets  into  a  force-free  shape  (ref.  8).  Note  that  the  form  of  the 
integrand  shown  in  the  body  integral  indicates  the  use  of  a  vorticity  distribution  to 
represent  a  thin  wing  in  some  proof- of-concept  calculations.  One  of  the  major 
generalizations  of  this  method,  to  be  initiated,  consists  of  replacing  this  body 
integral  with  the  UTSA  surface-panel  formulation  so  that  transonic  flow  over  bodies  of 
arbitary  shape,  including  vortex-type  separation,  can  be  calculated.  Other  planned 
improvements  include  (1)  replacing  the  vortex- 1 atti ce  model  used  in  the  wake  integral 
for  proof-of-concept  calculations  with  the  hybri d- vortex  formulation  (refs.  23  and 
24)  in  which  second-order  di stri buted- vorti ci ty  panels  are  used  to  compute  near-field 
influence,  reducing  to  zeroth-order  ( di screte- vortici ty )  elements  for  far-field 
influence,  (2)  shifting  the  linear  compressibility  term  xx  from  volume  integral 
to  surface  integral  by  solving 

v2<f>  -  M2<j>xx  =  Q  -  M2$xx  e  Qnonlin  (9) 

instead  of  equation  (7),  thereby  significantly  reducing  the  region  over  which  the 
volume  integral  needs  to  be  evaluated,  (3)  replacing  constant  source  strength  with 
linearly  varying  source  strength  in  the  volume  elements  and  introduing  a  threshold 
cutoff  value  for  the  integrand  of  the  volume  integral  to  terminate  integration  when  the 
integrand  diminishes  to  negligible  magnitude,  (4)  accelerating  convergence  of  the 
solution  by  possible  use  of  shock  fitting  (ref.  8),  (5)  accounting  for  entropy  changes 
across  shockwaves  (see,  e.g.,  ref.  12).  Code  development  for  unsteady  flow  is  in 
progress.  Research  on  suitable  configuration  of  these  codes  for  efficient  use  in 
computer-aided  interdisciplinary  design  will  be  a  continuing  activity. 

Completion  of  the  improvements  listed  above  should  provide  a  powerful  tool  for 
calculating  transonic  and/or  free-vortex  flows  around  arbitrary  aircraft  configurations 
with  sharp  leading  edges  or  with  specified  separation  line  locations.  Establishing  the 
separation  line  on  a  vibrating  wing,  however,  is  a  tough  viscous-flow  problem,  but  may 
be  amenable  to  treatment  by  the  seal ar- vector  potential  method  to  be  discusssed  below. 
The  importance  of  expediting  this  activity  should  be  underscored.  The  ability  to 
calculate  accurately  the  complicated  transonic  vortical  flows  around  highly  swept  wings 

and  complete  aircraft  at  high  angles  of  attack  is  a  key  problem  for  the  future 
development  of  highly  maneuverable  fighter  aircraft  and  is  already  needed  to  improve 
the  assessment  and  understanding  of  steady  and  transient  flight  loads  and  flutter 
problems  of  current  combat  aircraft.  It  should  be  especially  noted  that  vortex-type 
flow  separations  produce  typically  detrimental  effects  on  structural  loads  and  flutter. 

Figure  9  shows  the  calculated  velocity  field  and  shape  of  the  free-vortex  surface  in  a 
crossflow  plane  slightly  downstream  of  the  trailing  edge  of  a  delta  wing  with  vortex 
sheets  representing  thin  wakes  emanating  from  leading  and  trailing  edges  as  in  figure 
8.  The  volume  integral  (equation  (8))  has  not  been  included  for  this  incompressible- 
flow  calculation.  The  results  compare  quite  favorably  with  the  1 ow-Mach- number 
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experiments  of  Hummel  (ref.  25)  even  though  relatively  few  vortex  elements  were  used  in 
this  exploratory  calculation.  The  leading-edge  vortex  core  Is  clearly  defined  as  is 
the  incipient  deformation  of  the  trai 1 i ng-edge  vortex  sheet  into  a  trai 1 1 ng-edge  core 
with  rotation  opposite  to  that  of  the  leading-edge  core.  The  corresponding  spanwise 
distributions  of  lifting  pressure  ACp  are  shown  In  figure  10  for  crossflow  planes  at 
0.7  and  0.9  of  the  root  chord  aft  of  the  wing  apex.  Agreement  with  measured  values  Is 
very  good. 

Inclusion  of  the  volume  integral  (equation  (8))  permits  calculation  of  transonic  flow. 
Figure  11  shows  the  spanwise  distribution  of  upper- surf ace  pressure  Cpu  and  the  flow 
field,  including  a  captured  shock,  in  a  crossflow  plane  at  0.8  of  the  root  chord  aft  of 
the  apex  of  a  delta  wing  (ref.  8).  In  this  exploratory  calculation  the  vortex  sheet 
was  not  allowed  to  roll  up  enough  to  exert  Its  full  inductive  effect  on  the  wing 
surface  before  the  vorticity  was  transferred  Into  the  vortex  core.  If  an  additional 
quarter  turn  of  rollup  were  allowed,  the  pressure  peak  would  be  slightly  higher  and  a 
little  farther  outboard,  resulting  in  even  better  agreement  with  experiment.  In 
contrast,  the  pressure  peak  from  the  Euler  solution  is  considerably  weaker  and  farther 
outboard  than  the  experimental  peak  because  of  spatial  and  numerical  diffusion  in  the 
Euler  calculation. 

Structural  design  loads  do  not  occur  at  small - perturbati  on  conditions  but  at  limit 
load-factor  conditions  such  as  high  angle  of  attack.  Aeroelastic  deformations  are 
Important.  Wind-tunnel  results  may  be  of  questionable  accuracy  because  of  large  wall 
effects.  The  Important  influence  of  large  perturbation  conditions  and  free-vortex 
flows  on  structural  design  loads  is  typically  detrimental,  as  is  illustrated  by  the 
calculations  shown  in  figure  12  (from  ref.  26).  Even  if  the  linear  and  nonlinear 
spanwise  load  distributions  shown  were  compared  on  the  basis  of  same  total  normal  force 
(same  area  under  the  curves),  it  is  evident  that  the  effect  of  the  wing-tip  vortex  Is 
to  shift  the  load  outboard  and  hence  increase  wing  bending  movements. 

Linearized  aerodynamic  theory  Indicates  that  there  should  be  no  effect  of  angle  of 
attack  on  flutter  dynamic  pressure.  However,  a  detrimental  effect  typically  does  occur 
with  increasing  angle  of  attack  (see,  e.g.,  refs.  27  and  28).  If  adequate  flutter 
margins  are  to  be  maintained  when  angle  of  attack  is  not  near  zero,  the  degradation 
must  be  predictable.  Wind-tunnel  testing  of  sti ff ness- seal ed  flutter  models  is  not  the 
answer  because  they  are  typically  too  weak  to  sustain  more  than  very  small  static 
loads.  Figure  13  shows  experimental  variation  of  flutter  dynamic  pressure  with  angle 
of  attack  for  a  stiff  wing  that  was  spring  supported  (ref.  29).  The  initial  decline  in 
flutter  dynamic  pressure  between  0  and  7  deg  is  attributed  to  the  effect  of  the  tip 
vortex.  Confirming  calculations  by  methods  just  described  are  in  early  stages.  The 
drastic  decline  beyond  7  deg  Is  probably  caused  by  flow  separation  progressing  forward 
from  the  trailing  edge.  Prediction  of  that  behavior  will  require  solutions  of 
Navi er-Stokes  equations  as  discussed  below. 

Summarizing  the  status  of  1 ntegral -equati on  methods  for  vortex-type  (thin  wake)  flow 
separation:  The  hybri d- vortex  method  for  1 ow-Mach- number  steady  flow  (refs.  23  and  24) 
is  complete.  Computations  based  on  equation  (8)  for  steady  transonic  flow  with 
vortex-type  separation  and  shockwaves  have  been  demonstrated  (ref.  8),  and  the 
corresponding  unsteady  code  development  is  in  progress.  Major  generalizations  and 
improvements  in  efficiency  are  underway.  Further  developments  for  transonic  flow,  with 
or  without  vortex-type  flow  separation,  will  be  based  on  equation  (8). 

Scalar/Vector-Potential  Method:  When  viscous  influences  (other  than  thin  wakes  from 
TTTtT ng- surface  edges )  are  i mportant  --  for  example,  boundary- 1 ayer  effects  on 
control-surface  forces,  shock/boundary-layer  interaction,  or  flow  separation  from 
surfaces  (fig.  1)  --  solution  of  Navi er-Stokes  equations  in  some  form  is  required.  The 
approach  taken  here  1$  a  scalar/vector-potential  (SYP)  decomposition  of  the  velocity 
field  by  use  of  the  classical  Helmholtz  representation  of  a  vector  field  as  the  sum  of 
an  irrotational  part  and  a  solenoidal  part  (refs.  13  to  15).  Thus 

v  =  grad  <j>  +  curl  A  (10) 

where  $  is  again  the  scalar  potential  which  is  evaluated  by  the  methods  already 
described  herein,  and  the  vector  potential  A  is  related  to  the  vorticity  u  by 

v2A  =  -  u  =  -  curl  v  ( 11 ) 

The  vorticity,  In  turn,  is  governed  by  the  vorti ci ty-dynami cs  equation 

8i  (?)  -  ?  •  grad  *  =  }  cur1  A 

=  £  grad  T  x  grad  S  +  |  curl  ^  dlv  (T  +  pi) 


(12) 
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which  is  obtained  by  taking  the  curl  of  Navi er-Stokes  equation  for  general, 
three-dimensional,  unsteady,  compressible,  viscous,  heat-conducting  flow  (ref.  13)  to 
which  the  present  formulation  is  fully  equivalent.  In  equation  (12),  T  is  temperature, 
S  is  entropy,  and  T  is  stress  tensor. 

The  formulation  in  equations  (10)  to  (12)  appears  to  be  a  computationally  attractive 
alternative  to  direct  solution  of  the  Navi er-Stokes  equations  in  primitive  variables. 
Methods  of  this  type  have  been  used  for  a  long  time  for  viscous  incompressible  flow, 
but  they  have  not  proved  to  be  readily  general  1 zabl e  to  compressible  flow.  The  present 
formulation  is  quite  general  and  is  directly  applicable  to  compressible  flow.  Since 
the  outer  region  of  the  flow  about  an  aircraft  is  essentially  1 rrotati onal ,  an 
i ntegral -equati on  implementation  appears  to  be  an  especially  attractive  method  of 
solution.  The  initial  proof-of-concept  code  for  two-dimensional  incompressible  flow 
has  been  used  to  calculate  boundary  layers  on  a  flat  plate  (fig.  14),  flow  over  an 
airfoil,  and  separated  flow  around  a  rectangle  (fig.  15)  --  all  with  good  results 
(refs.  14  and  15).  In  particular,  the  calculation  of  flow  around  a  rectangle  (fig.  15) 
demonstrates  the  ability  to  calculate  flows  involving  large  regions  of  separation.  The 
resulting  velocity  field  shown  in  the  figure  illustrates  that  flow  separation  is 
predicted  very  satisfactorily  even  though  separation  was  not  imposed  by  any  artificial 
means  within  the  algorithm.  These  results  are  in  excellent  agreement  with  the 
finite-difference  results  of  ref.  31.  Application  to  a  circulation-control  airfoil  is 
being  initiated.  For  applications  to  turbulent  flows  this  method,  of  course,  requires 
a  good  turbulence  model  just  as  any  other  method  does.  In  addition  to  its 
computational  use,  the  SVP  formulation  has  also  generated  considerable  Insight  into  the 
relations  between  surface  boundary  conditions,  viscosity,  vorticity  and  its  diffusion 
(refs.  14  and  15). 

Current  activities  are  extending  the  proof-of-concept  code  to  three  dimensions  and  to 
compressible  flow.  The  types  of  applications  planned  Include  viscous  flow  over  lifting 
surfaces  with  and  without  control-surface  deflection,  lifting  surfaces  with  flow 
separation  from  edges  in  compressible  flow,  and  lifting  surfaces  with  separated  flow 
following  a  step  change  In  angle  of  attack. 

Summary  of  Integral-Equation  Activities 

The  activities  described  here  and  the  computational  capabilities  summarized  in  table  I 
indicate  that  completion  of  this  work  will  provide  efficient  and  unified  treatment  of 
flow  over  vehicles  having  arbitrary  shapes,  motions,  and  deformati ons  at  subsonic, 
transonic,  and  supersonic  speeds  up  to  high  angles  of  attack.  Moreover,  the 
computational  forms  of  the  equations  and  the  computational  capabilities  that  are 
emerging  appear  to  be  well  suited  for  repetitive  use  in  design  applications  as  well  as 
for  stand-alone  use.  As  pointed  out  previously,  the  UTSA  surface-panel  program  for 
attached  flow  may  contain  both  subsonic  and  supersonic  modules  in  a  single  program. 

Flow  complexities,  such  as  transonic  nonl 1 nearl ties ,  thin  wakes,  or  viscous  Influences, 
are  addressed  only  if  and  where  they  occur.  Thus,  if  the  volume-integral  module  Is 
Included  with  UTSA,  the  program  Implements  the  ful  1 -potenti al  equation  for  transonic 
nonlinear  attached  flow.  With  modification  for  shock-generated  entropy  change,  the 
program  can  apply  also  to  flows  with  shocks  of  finite  strength,  including  supersonic 
Mach  numbers  above  the  linear  range,  as  long  as  shock-generated  vorticity  Is  of  minor 
Importance.  If  the  hybrid- vortex  module  representing  the  free  vortex  sheets  is  also 
Included,  the  code  treats  transonic  flow  with  vortex-type  separation.  Finally, 
combination  of  the  vector  potential  with  these  scalar-potential  methods  (SVP 
formulation)  permits  the  formal  equivalent  of  Navi er-Stokes  solutions  for  high  angles 
of  attack  where  flow  separation  from  surfaces  may  occur  (for  example,  on  advanced 
fighter  aircraft  in  combat  maneuvers  and  in  highly  transient  conditions)  and  also  even 
for  low  angles  of  attack  when  control-surface  deflections  or  deflection  rates  are  large 
enough  or  shock  waves  are  strong  enough  to  cause  significant  boundary-layer  thickening 
or  separation.  The  latter  conditions  are  particularly  Important  for  generating  control 
forces  and  for  design  of  active  control  systems. 


FINITE-DIFFERENCE  METHODS 


The  goal  of  this  activity  is  to  develop  finite  difference  methods  that  can  be  used  for 
aeroelastic  analysis  of  complete  aircraft.  At  the  Langley  Research  Center,  efforts 
based  on  transonic  small  perturbation  (TSP)  potential  theory,  full  potential  theory 
and  the  Eul er/Navi er-Stokes  equations  are  underway. 

Transonic  Smal 1 -Perturbati on  Equation 


At  the  TSP  level,  development  has  progressed  on  two  fronts--(l)  extending  the 
capability  of  the  XTRAN3S  code  (ref.  32)  through  extensive  modification  and  (2) 
developing  a  new  program  (ref.  33).  Fig.  16  shows  a  sample  wing/fuselage  calculation 
made  possible  by  modifying  XTRAN3S  (ref.  34).  The  wing  has  an  RAE  101  airfoil  section, 
37  deg  leading  edge  sweep  angle,  aspect  ratio  (AR)  of  6,  and  taper  ratio  of  one-third. 
The  fuselage  is  a  sting-mounted,  axisymmetrlc  body  of  revolution  with  fineness  ratio 
l 1 engtn/maximum  diameter)  of  7.66.  Calculations  were  made  for  a  flexible  wing  at  free 
stream  Mach  number  (M)  of  0.91  and  mean  angle  of  attack  (a)  of  1  deg.  Fig.  16  (a) 
shows  the  wing  tip  deflection  as  a  function  of  time,  and  fig.  16  (b)  shows  the  griddinq 
used  to  represent  the  wing/fuselage.  Figs.  16  (c)  and  16  (d)  show  the  instantaneous 
Mach  number  contours  at  the  maximum  and  minimum  tip  deflections,  respectively.  The 
wing  motion  is  the  first  bending  mode  at  a  reduced  frequency  (k)  of  0.25.  The  contours 
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near  the  leading  and  trailing  edges  indicate  local  Mach  numbers  less  than  0.85,  and 
over  the  wing  chord  and  fuselage,  the  contours  Indicate  Mach  numbers  greater  than  0.95. 

A  new  TSP  code,  CAP-TSD  (Computational  Aeroel asti ci ty  P/ogram-l>ansoni  c  Small 
Disturbance)  (ref.  33)  has  been  developed  at  Langley.  It  solves  the  three-dimensional 
T3-D)  TSP  equation  using  an  approximate  factorization  (AF)  algorithm.  The  code  is 
significantly  more  efficient  than  methods  that  use  an  al ternati ng- di recti  on- i mpl i ci t 
(ADI)  solution  algorithm  and  can  be  used  for  aeroelastic  analysis  of  complete 
aircraft.  Fig.  17  shows  comparisons  of  CAP-TSD  (AF)  and  XTRAN3S  (ADI)  calculations 
with  experimental  data  for  a  rigid  F-5  wing  pitching  about  zero  mean  angle 
at  k  =  0.137  (ref.  35).  Pressures  are  shown  at  the  51  percent  span  station  (n). 

Upper  surface  pressures  are  shown  in  fig.  17  (a),  and  lower  surface  pressures  are 
shown  in  fig.  17  (b).  Both  sets  of  calculations  show  good  agreement  with  the  measured 
data  and  are  nearly  identical  to  each  other.  The  primary  difference  is  that  the 
AF  solution  requires  only  ten  percent  of  the  computer  resources  used  in  the  ADI 
calculation.  CAP-TSD  has  been  used  to  compute  steady  flow  past  a  wl ng/f usel age/tall 
model  that  was  tested  at  DFVLR .  In  the  test,  M  =  0.2,  a  =  0.15  deg.  The  wing  is 
rectangular  with  RAE  101  airfoil  sections  (9  percent  thickness  ratio)  and  a  full-span 
aspect  ratio  of  6.  The  horizontal  tall  is  rectangular  with  RAE  101  airfoil  sections 
(12  percent  thi cknessrati o )  and  full-span  aspect  ratio  of  3.  The  fuselage  Is  an 
axisymmetric  body  of  revolution  with  fineness  ratio  of  9.75.  The 
mathematical  representation  of  the  model  is  shown  in  fig.  18  (a).  Comparisons  of 
computed  and  measured  pressures  on  the  wing  and  tail  are  shown  in  fig.  18  (b).  Fig.  18 
(c)  shows  comparisons  of  computed  and  measured  pressures  on  the  fuselage.  In  all 
cases,  the  computations  and  experiments  are  in  good  agreement. 

Potential  flow  theory  has  been  shown  to  give  highly  erroneous  and  even  multivalued 
results  when  shock  waves  are  in  the  flow  field  (refs.  36,  37).  This  is  because 
shock-generated  entropy  is  not  modeled  in  potential  flow  formulations.  A  method  for 
modeling  noni sentropi c  effects  in  2-D  TSP  theory  was  developed  by  Fuglsang  and  Williams 
(ref.  12)  and  extended  to  three  dimensions  by  Gibbons  et  al .  (ref.  38).  The 
noni sentropi c  formulation  was  implemented  by  modifying  the  streamwise  flux  in  the  TSP 
equation  to  account  for  entropy  jumps  across  shock  waves.  This  alleviates  the 
phenomena  of  multiple  solutions  and  highly  Inaccurate  loading  predicted  by  Isentropic 
potential  methods.  Fig.  19  shows  an  example  of  calculated  lift  as  a  function  of  angle 
of  attack  for  isentropic  and  noni sentropic  formulations.  An  Euler  calculation  is 
Included  at  one  deg  angle  of  attack.  Without  the  noni sentropi c  corrections,  the 
calculated  lift  is  too  large.  When  the  corrections  are  included,  the  calculated  lift 
is  less  and  agrees  with  the  Euler  calculation  at  the  point  where  such  data  is 
avai 1 abl e. 

Ful 1 -Potenti al  Equation 

A  method  for  modeling  shock-generated  entropy  in  the  unsteady  full-potential  (FP) 
formulation  has  been  developed  at  Langley  (ref.  39).  The  method  is  an  extension  of  the 
steady-flow  method  of  Hafez  and  Lovell  (ref.  40).  Fig.  20  shows  the  Instantaneous 
pressures  on  an  NACA  0012  airfoil  oscillating  in  pitch  about  its  quarter  chord.  In 
this  case,  M  =  0.755,  a(t)  =  (0.016  +  2.51s1n(kt))  deg,  and  k  (based  on  semichord)  = 
0.0814.  Fig.  20  (a)  shows  a  comparison  of  the  isentropic  full -potenti al  method,  a  TSP 
method  (denoted  by  "TSD"  on  the  figure)  (ref.  41),  and  experimental  data  (ref.  42). 

Fig.  20  (b)  shows  Isentropic  and  noni sentropi c  full  potential  methods,  along  with  the 
measured  data.  When  entropy  corrections  are  used,  more  accurate  modeling  of  the  shock 
is  obtained,  and  agreement  with  the  measured  data  Is  improved. 

Currently,  efforts  are  underway  to  extend  the  full  potential  method  to  three 
dimensions.  The  proposed  method  will  have  the  capability  to  do  aeroelastic  analysis 
for  flows  in  the  subsonic,  transonic,  and  supersonic  speed  ranges. 

Eul er/Navi er-Stokes  Equations 

Research  is  being  conducted  to  develop  methods  for  obtaining  time-accurate  unsteady 
solutions  of  the  Eul er/Navi er-Stokes  equations.  These  methods  are  used  to  solve  the 
Navi er-Stokes  equations,  with  solutions  of  the  Euler  equations  obtained  by  turning  off 
the  viscous  terms.  The  existing  codes  are  used  to  march  In  time  to  steady  state  and 
currently  are  being  modified  to  be  made  time-accurate.  An  effort  is  underway  to 
correlate  Navi er-Stokes  calculations  with  pressure  data  measured  in  the  Langley 
Research  Center  0.3-meter  Transonic  Cryogenic  Tunnel.  Unsteady  transonic  pressures 
were  measured  on  a  14- percent- thi ck  supercri ti cal  airfoil  at  cryogenic  temperatures  for 
M  between  0.65  and  0.74  and  at  Reynolds  numbers  based  on  airfoil  chord  (Rc)  between  6 
million  and  35  million  (ref.  43).  The  model  is  shown  in  fig.  21  (a).  The  open 
symbols  in  fig.  21  (b)  show  test  conditions  where  the  effects  of  frequency  on  the 
unsteady  pressures  were  examined,  and  the  solid  symbols  show  where  the  effects  of 
frequency  and  amplitude  were  studied.  Steady  pressure  distributions  at  Rc  =  6  million 
(open  symbols)  and  at  Rc  =  30  million  (solid  symbols)  are  shown  in  fig.  21  (c).  The 
results  of  this  study  may  be  used  to  determine  when  various  forms  of  the  Navi er-Stokes 
equations  can  be  used  (e.  g.,  thin-layer  or  the  full  equations)  and  may  be  used  to 
evaluate  methods  that  couple  viscous  flow  models  with  inviscid  methods. 

An  implicit  upwind  code  that  can  be  used  to  calculate  massively  separated  2-D  flows 
(ref.  44)  is  available  for  use  in  the  correlation  study.  It  uses  van  Leer  flux-vector 
splitting  and  is  first-order  accurate  in  time  and  second-order  accurate  in  space.  The 
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method  can  be  used  for  time  marching  to  steady-state  solutions,  but  It  Is  not 
time-accurate  and  cannot  model  unsteady  flow.  An  example  of  the  capability  of  this 
code  is  shown  In  fig.  22.  Even  though  the  code  Is  not  time  accurate,  the  example  shows 
the  variation  in  loads  that  can  occur  when  marching  to  steady  state.  It  shows 
computation  of  laminar  flow  about  a  12-percent- thick  Joukowski  airfoil  at  53  deg  angle 
of  attack  at  minimum  lift  (fig.  22  (a)).  Increasing  lift  (fig.  22  ( b ) ) ,  and  maximum 
lift  (fig.  22  ( d ) )  for  a  Strouhal  number  of  0.166.  The  time  history  of  the  lift  is 
shown  In  fig.  22  (c).  At  minimum  lift,  a  strong  counterclockwise  vortex  has  just  been 
shed  from  the  trailing  edge.  There  is  a  weak  region  of  clockwise  vortlcity  one-half 
chord  behind  the  airfoil  and  a  strong  region  of  clockwise  vortlcity,  which  persists 
throughout  the  entire  cycle,  at  the  leading  edge.  The  flow  near  the  upper  surface  is 
broken  into  several  cells  of  vorticity.  As  the  lift  Increases,  the  trailing  edge 
vortex  breaks  away  and  weakens  as  the  region  of  clockwise  vortlcity  behind  the  airfoil 
increases  in  strength.  The  vorticity  near  the  upper  surface  becomes  predominantly 
counterclockwise.  Current  efforts  are  aimed  at  extending  this  capability  to 
time-accurate  analysis  of  2D  unsteady  flows. 

A  steady-flow  Navi er-Stokes  method  (ref.  45)  Is  available  for  analysis  of  3-D  wings. 

The  convective  and  pressure  terms  are  upwind  differenced,  using  a  flux-vector  splitting 
method.  The  shear  stress  and  heat  transfer  terms  are  centrally  differenced.  The 
resulting  algorithm  is  second-order  accurate  in  space.  An  implicit,  spatially  factored 
algorithm,  which  is  fully  vectorized  for  the  Control  Data  Corporation  VPS  32,  is  used 
to  provide  efficient  solutions.  Efforts  also  are  underway  to  extend  this  capability  to 
3D  time-accurate  analysis.  An  example  of  calculations  made  using  this  method  is  shown 
in  fig.  23.  It  shows  the  contours  of  the  calculated  total  pressures  on  an  AR  =  1  delta 
wing  (75  deg  leading-edge  sweep  angle)  at  M  =  0.3  and  a  -  20.5  deg.  Solutions  were 
obtained  by  marching  in  time  to  steady  state.  Primary  and  secondary  vorticies  are 
evident  on  the  upper  surface.  Close  examination  of  the  surface  velocities  indicates  a 
tertiary  separation  outboard  of  the  secondary  vortices.  This  type  of  flow  field  also 
was  observed  in  a  related  experiment  (ref.  25). 

Fig.  24  shows  contours  of  measured  total  pressures  on  the  same  delta  wing  at  the  same 
flow  conditions  used  in  the  previously  mentioned  Navi er-Stokes  calculations.  These 
steady-flow  data  were  measured  in  the  Basic  Aerodynamic  Research  Tunnel  at  Langley 
(ref.  46)  and  show  good  agreement  with  the  calculated  data  in  fig.  23. 

Summary  of  Finite-Difference  Activities 

Areas  of  current  activity  in  the  development  of  f i ni te- di f f erence  methods  are  shown  in 
fig.  25.  Methods  based  on  TSP  theory  are  being  applied  to  configurations  as  complex  as 
complete  aircraft.  Unsteady  ful 1 -potenti al  methods  are  being  developed  for  3-D 
configurations  with  the  goal  being  to  use  such  methods  for  aeroelastic  analysis  of 
complete  aircraft.  Steady,  time-marching  Navi er-Stokes  methods  are  available  for 
airfoil  and  for  isolated  wings.  Those  methods  are  currently  being  made  time-accurate. 


CONCLUDING  REMARKS 

Some  problems,  progress,  and  plans  in  the  development  of  steady  and  unsteady 
computational  aerodynamics  for  use  in  aeroelastic  analysis  and  design  have  been 
reviewed.  The  primary  focus  has  been  on  applications  to  (1)  vehicles  having  arbitrary 
shapes,  motions,  and  deformations,  (2)  appropriate  design  and  operating  conditions, 
especially  for  transonic  speeds  and  high  angles  of  attack,  (3)  efficient  computation  of 
aerodynamics  and  aeroelastic  behavior  for  both  design  and  analysis.  Current  and  future 
activities  have  been  highlighted. 
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Table  1. -Summary  of  Integral-Equation  Activities. 


— - - M  Range 

a  Range  — — 

Subsonic 

Transonic 

Supersonic 

Low 

(attached  flow) 

w/large  control  deflection 

UTSA 

SVP 

Nonlinear  UTSA 

SVP 

UTSA 

Moderate 

(vortex  separation) 

w/large  control  deflection 

UTSA  + 
Hybrid  vortex 

SVP 

Nonlinear 
UTSA  -f 
Hybrid  vortex 

SVP 

Large 

(separated  flow) 

w/  or  w/o  control  deflection 

SVP 

SVP 

SVP 

SVP 

1980 


1990 


Sub/supersonic  integral  eq.  method  for  complete  a/c 
Linear  Nonlinear 
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?EF-an&  3D  Transonic  methods 

Mn\/i  or-QtnLoc  on  i 

x‘a 

INUVlcl  o LUIsco  t.  L| »  / 

;  Viscous  flow  effects  £ 

Edge  separation;  B.  L.  effect;  shock/B.  L.; 

Surface  sep. 

_ \ 

Deleted  from  original  plan  I  I  Added  to  original  plan 


Goal :  Validated  computdtional  methods  for  evaluating  steady  dnd 
unsteady  loads  on  aircraft  having  arbitrary  shapes,  motions,  and 
deformdtions  (including  control  surfaces)  in  subsonic,  transonic, 
and  supersonic  flow  up  to  high  angles  of  attack 

Justification:  These  methods  dre  needed  to  evoluate  ond  study 
structurdl  loads,  aerodynamic  coefficients,  stability  character¬ 
istics,  dynamic  loads,  and  flutter  in  the  analysis  and  design  of 
advanced  high-performance  aircraft,  including  fighters  capable 
of  supermanuverability 


Fig.  1  -  Integral -equation  program  in  unsteady  aerodynamics 
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Fig.  2  -  Shuttle  orbiter  flutter  analysis 
by  integral -equation  method 


Fig.  3  -  Unsteady  surface  pressures  on 
clipped-delta  wing  at  Mach  number  0.4, 
reduced  frequency  0.66 


Fig.  4  -  Effect  of  angle  of  attack  on 
steady  pressures  for  clipped-delta 
wing  with  6-percent  biconvex  airfoil 
at  Mach  number  0.4,  Reynolds  number  2.xl0^ 


Fig.  5  -  Steady  pressures  near  root  of 
aspect-ratio  6  rectangular  wing  with 
NACA  0012  airfoil  at  Mach  number  0.82 

a-0 


1-14 


0 


o  Integral  equation 

o  Experiment 
r  Sonic 

pressure  ro$8o 

-  o 


o 

0 


6> 


®  Fraction  of  chord  ^ 


AC 


Fig.  6  -  Steady  pressures  for  aspect- 
ratio  4  rectangular  wing  with  6%  bi¬ 
convex  airfoil  at  Mach  number  0.908, 

a=0 


Fig.  7  -  Unsteady  pressures  calculated 
with  SUSAN  code  for  aspect- 
ratio  5  rectangular  wing  with  NACA 
64A006  airfoil  pitching  at  reduced 
frequency  0.06,  Mach  number  0.875 
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Fig.  8  -  Sketch  of  vortex  sheets  separating 
from  wing  edges 


Fig.  9  -  Flow  field  behind  aspect- 
ratio  1  delta  wing  at  <*=20.5° 
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Fig.  10  -  Steady  pressures  on  aspect- 
ratio  1  delta  wing  at  <*=20.5° 


Fig.  11  -  Steady  pressures  on  aspect- 

ratio  1.5  delta  wing  at  Hach  nucber  0.7, 
0*15.°,  x/cr=0.8 


Dynamic 
press. ,psf 


l - i - 1 - 1 _ l _ i _ _ _ ( 

0  2  4  6  8  10  12 

a,  deg 


Fig.  12  -  Steady  spanwise  variation  of  local  Fig.  13  -  Effect  of  angle  of  attack  on 
nornal-force  coefficient  for  aspect-  flutter  of  aspect-ratio-6  rectangular 

ratio  1  rectangular  wing  at  a=19.4°  wing  with  NACA  64A010  airfoil 
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Fig.  15  -  Velocity  field  around  a  rectangle  in  incompressible  flow 
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(a)  Tip  deflection 


(b)  Wing/fuselage  gridding 


(c)  Maximum  tip  deflection  (A) 


(d)  Minimum  tip  deflection  (B) 


Fig.  16  -  Unsteady  transonic  calculations  for  RAE  wing/fuselage  at  Mach 
number  0.91,  a  =  1  deg 


Fig.  17  -  First  harmonic  of  unsteady  pressures  on  F-5  wing  at  Mach  number 
0.9,  a  =  0  deg,  k  =  0.137 
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(a)  Mathematical  model 


Fig.  19  -  Calculated  lift  on  a  rectangular 
wing  with  NACA  0012  sections  at  Mach 
number  0.82 
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(b)  Wing  and  tail  pressures 
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(a)  Isentropic  vs.  experiment 


(c)  Fuselage  pressures 

Fig.  18  -  Steady  flow  computations  for  DFYLR 
model  at  Mach  number  0.2,  a  =  0.15  deg 


Xfr. 


(b)  Isentropic  and  noni sentropi c  vs. 
experiment 


Fig.  20  -  Unsteady  pressures  on  an  NACA  0012 
airfoil  at  Mach  number  0.755,  a(t)  = 

(0.016  +  2.51sin(kt))  deg,  kt  =  168  deg 
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(c)  Steady  pressure  distributions  at  Mach 
number  0.72,  a  =  1.5  deg 


Fig.  21  -  Study  of  Reynolds  number  effects  on  unsteady  pressures 


(c)  Lift  history  (d)  Maximum  lift 


Fig.  22  -  Computation  of  two-dimensional  flow  separation 
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Fig.  23  -  Calculated  total  pressure  contours 
on  a  75  deg  delta  wing  at  Mach  number 
0.3,  a  =  20.5  deg,  Reynolds  number  =  0.95 
■ill  ion 


Fig.  24  -  Measured  total  pressure  contours 
on  a  75  deg  delta  wing  at  Mach  number  0.3, 
a  =  20.5  deg,  Reynolds  number  =  1  million 
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Fig.  25  -  Summary  of  finite-difference  activities 
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SUMMARY 

Unsteady  airloads  in  the  transonic  speed  range  for  the  application  in  flutter  calculations  can  only  be 
predicted  by  using  nonlinear  flow  equation  codes.  In  this  paper  a  typical  section  with  the  MBB-A3  profile 
is  used  as  theoretical  model  for  the  flutter  investigations.  Flutter  calculations  are  presented  using  the 
classical  method  in  the  frequency  domain  as  well  as  a  structural -aerodynamic  coupling  procedure  in  the  time 
domain.  The  unsteady  aerodynamic  coefficients  which  are  used  for  these  flutter  investigations  are 
calculated  for  several  M-numbers  and  incidences  with  the  TSP-equation  without  and  with  quasiunsteady 
boundary  layer  and  with  the  time  linearized  TSP  equation  without  and  with  steady  boundary  layer  and  for 
comparison  reasons  with  the  Doublet  Lattice  method. 

The  boundary  layer  influences  essentially  the  unsteady  aerodynamic  pressure  distributions  and 
coefficients.  Mainly  for  small  reduced  frequencies,  the  coefficients  calculated  with  TSP  deviate  strongly 
from  the  coefficients  calculated  with  Doublet  Lattice,  but  the  quasiunsteady  boundary  layer  reduces  these 
differences  significantly.  The  results  achieved  with  linearized  TSP  behave  similarly,  but  with  smaller 
differences  compared  to  TSP.  The  flutter  speeds  calculated  with  TSP  airloads  show  the  expected  dip  in 
comparison  to  the  Doublet  Lattice  airloads.  The  dip  over  the  transonic  M-number  range  is  reduced  in  the 
range  of  higher  flutter  frequencies.  But  the  boundary  layer  also  reduces  the  dip  considerably.  The 
boundary  layer  also  can  raise  the  flutter  speed  above  the  Doublet-Lattice  values  for  M-numbers  higher  than 
the  design  M-number  when  the  boundary  layer  influence  changes  the  flutter  frequency. 

The  flutter  results  obtained  in  the  frequency  and  time  domain  differ  less  from  each  other.  Nevertheless 
the  time  marching  flutter  simulation  takes  the  nonlinearity  effects  of  the  airloads  dependent  on  the 
vibration  amplitude  essentially  better  into  account  than  the  classical  flutter  solution. 


1.  INTRODUCTION 

Since  the  publication  of  Farmer  and  Hanson  [1]  it  became  apparent  that  flutter  calculation  results  in 
the  transonic  speed  range  were  nonconservative  in  comparison  to  the  realistic  ones  when  obtained  with 
unsteady  airloads  calculated  by  linear  theories,  especially  for  wings  with  modern  profiles.  Since  that 
time  there  has  been  extensive  progress  in  transonic  unsteady  computational  aerodynamics  for  application  to 
flutter  analysis  12  ]  .  But  the  application  of  inviscid  transonic  codes  for  flutter  calculations  did  not 
improve  the  situation.  In  the  range  of  small  shocks  the  description  of  the  inviscid  unsteady  flow  by  an 
Euler  code  offered  no  advantages  over  the  TSP  code  (2  J  ,  (31.  More  realistic  flutter  results  were 
obtained  by  taking  into  account  viscous  effects  on  the  unsteady  airloads  by  profile  contour  changes  due  to 
boundary  layer  displacement  thickness  and  shock-boundary  layer  interaction.  In  unsteady  aerodynamics 
viscous  effects  are  at  least  as  important  as  in  steady  aerodynamics,  because  they  not  only  affect  the  shock 
position  and  strength,  but  also  the  extent  of  shock  movement. 

In  contrast  to  unsteady  pressure  distributions  of  linear  aerodynamic  theories  the  unsteady  pressure 
di stributions  obtai ned  in  the  transonic  speed  range  depend  on  the  mean  incidence  of  the  profile  and  of  the 
amplitude  of  its  motion.  Additionally  the  unsteady  pressure  distribution  induced  by  harmonic  motion  of  the 
profile  does  not  behave  harmonically  but  only  periodically,  i.e.  such  an  unsteady  pressure  distribution 
contains  higher  harmonics  for  a  harmonic  motion  of  the  profile. 

The  flutter  solutions  in  the  frequency  domain  employ  only  the  first  harmonic  of  the  airloads  for  each 
vibration  mode  in  an  eigenvalue  analysis  to  determine  the  flutter  point.  By  using  time  marching  codes  for 
producing  these  airloads  the  airfoil  is  forced  to  oscillate  in  the  prescribed  mode  at  the  desired  frequency 
until  the  transients  have  decayed  and  a  periodic  solution  is  achieved.  After  Fourier  analysis  of  the  time 
histories  only  the  first  harmonic  is  used  in  the  flutter  calculation.  The  results  for  different 
oscillation  modes  are  superposed  in  the  solution  process  to  represent  the  airloads  due  to  the  flutter  mode. 

In  the  time  domain  flutter  simulation  the  structural  and  aerodynamic  equations  are  coupled  and  the 
complete  system  is  time-marched  from  a  prescribed  excitation.  The  stability  is  assessed  from  the  growth  or 
the  decay  of  the  transient  response.  With  this  method  the  higher  harmonics  of  the  airloads  and  the 
influence  of  different  forms  of  disturbances  on  the  stability  behaviour  can  be  taken  into  account  and 
superposition  of  modes  is  no  longer  necessary. 

In  this  paper  flutter  calculation  results  are  shown  for  a  typical  section  with  the  MBB-A3  profile,  see 
also  14],  (7  ].  For  the  calculation  of  the  airloads  a  finite  difference  TSP  code  is  used  which  was 
developed  by  the  ONERA  [5],  For  the  flutter  calculation  in  the  frequency  domain  the  airloads  were 
calculated  for  various  reduced  frequencies  with  the  TSP-equation  as  well  as  with  a  time  linearized 
TSP-equation.  Viscous  effect  have  been  taken  into  account  by  an  integral  boundary  layer  method  [6], [8]. 
After  each  time  step  of  the  TSP  calculation  a  boundary  layer  calculation  was  performed. 
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2.  AERODYNAMIC  TOOLS 

For  the  inviscid  flow  computation  the  TSP-equation  for  the  velocity  potential: 


M2k2pu  +  2M2kpxt= 

n  <J>  2  ^  +  <$> 

2  ^  x  J  yy 

(i) 

with 

X  =  (y  +  1 )  M2 

(la) 

is  solved  using 

an  ADI  Finite-Difference-Scheme.  The  procedure  was 

establ  i  shed  by  [  5  ]. 

The  boundary 

conditions  on  the  airfoil:  y=f(x,t)  and  on  the  wake 

are: 

=  fx  +  ft  + 

on  the  airfoil 

(2) 

A  <t>y  =  A  dx 

(3a) 

A(*x  +  4.t)  =  0 

across  the  wake 

(3b) 

d  means  the  boundary  layer  displacement  thickness  and  A  indicates  the  jump  (+y-*0)  -  ( — y  -*  0 ) 


For  the  pressure  derivative  the  linearized  expression 

C p  =  — 2  (  4>x  +  <t>,)  (4) 

is  used.  Across  the  wake 

ACp  =  0  (4a) 

For  the  time  linearized  version  of  the  TSP-equation: 

M2k2^tt  +  2M2k«5)rt=  g-x(  +  X*x°  9x)  +?„  (5) 

with  the  steady  potential  <J>°  the  same  solution  procedure  is  applied  as  for  the  TSP-equation. 

A  steady  integral  method  [6]  is  used  for  the  boundary  layer  calculation.  Usually  for  the  calculation  of 
the  boundary  layer  displacement  thickness  an  iterative  invi scid-vi scid  coupling  procedure  is  used.  For 
this  procedure  the  velocity  distribution  of  the  last  time  step  is  used  as  boundary  condition  for  the 
calculation  of  the  next  time  step  of  the  displacement  thickness  (Fig.  1).  But  in  the  vicinity  of  the 
inviscid  separated  flow  region  one  of  the  boundary  layer  equations  become  singular.  Therefore  for  the 
separated  flow  region  an  indirect  solution  of  the  boundary  layer  equation  is  used;  that  means  an  estimated 
boundary  layer  displacement  thickness  is  calculated  by  an  inverse  procedure  and  in  each  iteration  step  the 
velocity  distribution  of  the  inviscid  region  is  compared  with  the  velocity  distribution  of  the  viscid 
region  by  an  relaxation  formula  (Fig.  2).  In  this  paper  the  Carter  relaxation  formula  is  applied. 

For  the  calculation  of  the  unsteady  pressure  distribution  including  the  boundary  layer  the  following 
three  phases  for  the  time  marching  process  are  used: 

1.  For  the  steady  pressure  distribution  calculation  with  the  TSP-equation  the  real  profile  is  blown  up  from 
the  flat  plate  in  about  80  iteration  steps. 

2.  After  that  the  steady  boundary  layer  displacement  thickness  is  added  in  an  alternating  process  to  the 

geometrical  profile  thickness.  This  process  needs  30  up  to  70  time  steps. 

3.  After  the  steady  phase  is  finished  the  harmonic  excitation  of  the  profile  is  started.  For  each  time 
step  only  one  vi scid-invi scid  calculation  is  made,  as  long  as  the  boundary  layer  calculation  does  not 
diverge.  No  unsteady  terms  are  applied  in  the  boundary  layer  equation. 

Fig.  4  shows  the  growth  of  the  boundary  layer  displacement  thickness  with  the  increasing  number  of 
iteration  steps  for  the  steady  case  of  the  MBB-A3  profile  (Fig.  3)  for  M=0.8,  a=-0.2°  and  Re=17.7  10s.  The 
ramp  in  the  boundary  layer  displacement  thickness  in  front  of  the  shock  position  is  visible.  The  triangle 
shows  the  laminar-turbulent  transition  point.  The  cross  indicates  the  chordwise  downstream  position  at 
which  the  direct  boundary  layer  calculation  was  replaced  by  the  indirect  one. 

Fig.  5  shows  the  steady  and  unsteady  pressure  distribution  of  a  TSP  calculation  without  and  with 
boundary  layer.  The  steady  parameters  are  the  same  as  in  Fig.  4  (  M=0.80,  a=-0.2°  Re=17.7  106  )  for  the 
reduced  frequency  of  k=0 . 1 04 .  The  result  of  an  Euler  calculation  without  boundary  layer  for  the  same  flow 
parameters  is  added.  The  results  for  the  steady  and  unsteady  pressure  distributions  for  TSP  without 
boundary  layer  and  Euler  agree  well.  The  large  influence  of  the  boundary  layer  on  the  pressure 
distribution  is  clearly  visible. 


Aeroelastic  Equation  of  Motion 

The  classical  flutter  calculation  is  based  on  a  modal  representation  of  the  structure.  The  elastic 
displacements  are  developed  in  a  series  of  structural  eigenmodes  : 

f  (  X,  t  )  =  t  %  (  X)  q,(t) 
i-i 


(6) 
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For  a  two-degree  of  freedom  case  this  representation  is  exact  for  N=2.  The  structural  equation  of  motion 
then  reads: 


Mi+  =  Z.(t) 


q 

M 

K 

Z 


generalized  coordinates 
generalized  masses 
generalized  stiffness 
generalized  aerodynamic  forces 

Z;(t)=  -LpV2  f  Ap(t)  +  i  dF 


(7) 


(8) 


For  the  flutter  calculation  in  the  frequency  domain  only  the  first  harmonic  of  the  airloads  induced  by  a 
pure  harmonic  motion  is  taken  into  account.  The  generalized  aerodynamic  forces  are  supposed  to  be  linear 
functions  of  the  generalized  coordinates  and  thus  can  be  superposed: 


Zi(t)=  E 
i- 1 

Ajj  <*>) 

(9) 

with 

Ag=  ipv\ 

/ ApjCw)^,  dF 

(9a) 

In  this  way  the  differential 

equation  (7)  can  be 

transformed  into  an  eigenvalue  problem 

[-x2m  +  k 

i 

> 

/'"'N 

e 

i— i 

.Ol 

II 

O 

(10) 

with  iX  =  <5  +\g)  (10a) 

The  implicit  dependence  A(u)  only  allows  an  iterative  solution  process.  In  this  paper  the  flutter 
calculations  are  performed  with  the  p-k  method  19]. 

For  the  flutter  simulation  in  the  time  domain  equation  (7)  is  solved  directly.  This  approach  is  not 
restricted  by  linear  assumptions  of  the  aerodynamic  forces.  In  this  paper  flutter  simulation  results  are 
obtained  using  a  Newmark  integration  [10  1.  The  equation  of  motion  (7)  then  reads: 


q  n+1  =  [  aQM  -4-  K  ]  [  Z  n+1  +  o0M  q  +  o^M  q n  4-  a3M  q  ]  t 

( 1 1  a ) 

qn+1  =  a0  (q n+1  -  q")-  a2qn  -  a3  q" 

(lib) 

•n+1  —  An  /■,  'An  . l  n+1 

q  -q  +a6q  +  a7q 

(11c) 

with 

a,-  -  integration  constants 

n  -  time  step  index 

The  aerodynamic  forces  at  the  time  level  n+1  are  still  unknown  when  solving  equation  (11a).  Thus  those 
of  time  level  n  are  introduced.  This  causes  the  phase  shift  of  one  time  step  which  can  be  neglected  for 
sufficient  small  time  steps. 

Frequency  and  damping  estimates  are  determined  fitting  the  time  responses  by  complex  exponential 
functions  of  the  form: 

q(x.t)  -  a0  +  2  e6it  [  cos  c^t  +  b-,  sin  w-t  ]  (]2) 

A  least  squares  curve  fitting  procedure  is  employed  following  an  optimization  algorithm  of  Jacob  (111. 


3.  RESULTS  AND  DISCUSSION 
Overview 


Flutter  calculation  were  established  for  a  typical  section  with  an  MBB-A3  profile  (Fig.  3).  Steady 
pressure  distributions  in  the  M-number  range  between  M=0.756  and  M=0.82  were  calculated  with  the 
TSP-equation  for  various  incidences.  Unsteady  pressure  distributions  for  the  same  M-number  and  incidence 
range  were  calculated  in  a  reduced  frequency  range  between  0.05  and  0.5  for  pitch  and  heave  motion.  It  was 
striking  that  the  solutions  for  unsteady  pressure  distributions  for  M-numbers  higher  than  M=0.8  diverged  in 
most  cases  although  the  steady  shock  was  not  positioned  at  the  trailing  edge  of  the  profile.  Solutions 
with  the  time  linearized  TSP-equation  for  selected  cases  were  also  established.  For  these  cases  the 
necessary  steady  potential  was  calculated  with  the  TSP-equation.  For  some  cases  the  steady  and  unsteady 
pressure  distributions  were  calculated  including  the  steady  and  quasi -unsteady  boundary  layer.  For  the 
calculation  of  the  unsteady  pressure  distribution  with  the  linearized  TSP-equation  including  the  boundary 
layer  only  the  boundary  layer  for  the  steady  potential  was  taken  into  account.  The  following  abbreviations 
are  used  for  the  aerodynamic  results. 


D.L . Doublet  Lattice 

TSP . Transonic  small  perturbation 

Tin.  TSP . Time  linearized  TSP 

vise.  TSP . TSP  with  quasi  unsteady  boundary  layer 

vise.  lin.  TSP . Time  linearized  TSP  with  steady  boundary  layer 
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Because  the  solutions  with  the  term  cj>^  in  the  TSP-equation  showed  disturbances  this  term  was  neglected. 
This  neglect  is  without  an  influence  on  the  flutter  results,  because  the  reduced  flutter  frequencies  of  the 
performed  flutter  calculations  were  smaller  than  k=0.3.  From  the  first  harmonic  of  the  unsteady  pressure 
distributions  the  force-  and  moment-coefficients  were  derived.  The  unsteady  aerodynamic  coefficients  with 
D.L.  and  for  the  plate  with  a  purely  elliptic  flow  field  with  TSP  were  calculated  for  reasons  of 
compari son. 


Steady  Aerodynamic  Results 


Fig.  6  shows  the  steady  pressure  distributions  for  M=0.765  anda=1.5°  calculated  with  the  TSP  equation 
including  the  boundary  layer  effect  in  comparison  with  A.R.A.  measurements.  The  calculated  pressure 
distribution  fits  the  measured  one  well.  The  pressure  distribution  calculated  by  TSP  (  matching  the 
M-number  and  the  incidence  )  does  not  fit  as  well  to  the  measured  one  as  the  vise.  TSP  solution  which  was 
not  matched.  But  the  boundary  layer  effect  on  the  steady  pressure  distribution  can  only  approximately  be 
represented  by  inviscid  calculations  and  changed  parameters. 

Fig.  7  shows  the  steady  pressure  distribution  for  M=0.78  and  a =1.3°  calculated  with  TSP  and  vise.  TSP. 
The  pressure  distribution  without  boundary  layer  has  a  more  downstream  shock  position  and  a  higher  pressure 
rise  at  the  shock  than  the  pressure  distribution  including  the  boundary  layer,  resulting  in  a  reduction  of 
steady  lift  in  the  latter  case.  The  vise.  TSP  solution  can  approximately  be  matched  by  TSP  when  reducing 
the  incidence  to  a  =0.75°. 


Unsteady  Aerodynamic  Results 


Fig.  8  and  Fig.  9  shows  the  unsteady  pressure  distribution  for  M=0.756  and  a=1.3°  with  a  pitch 
amplitude  of  0.5°  for  the  reduced  frequency  range  0.1  to  0.5  calculated  with  TSP  and  vise.  TSP 

respectively.  The  peaks  in  the  unsteady  pressure  distributions  on  the  upper  side  near  the  steady  shock 

positions  are  strongly  decreased  and  shifted  upstream  by  the  boundary  layer  effect.  But  even  the  levels  of 
the  pressure  fore  and  aft  of  the  peaks  are  changed. 

Fig.  10  shows  the  unsteady  pressure  distribution  calculated  for  M=0.78,  a-1.3°  and  a  pitch  angle  of 
0.5°  and  a  reduced  frequency  of  k=0.2  with  TSP  and  vise.  TSP.  The  unsteady  pressure  distribution  for 
a =0.75°  calculated  without  boundary  layer  is  shown  for  comparison.  The  latter  does  not  fit  the  unsteady 
pressure  distribution  calculated  with  TSP  including  boundary  layer  as  well  as  the  steady  pressure 
distribution  do  under  the  same  conditions  (see  also  Fig.  7). 

Fig.  11  shows  the  unsteady  pressure  distribution  for  M=0.78;  a  =0.5°;  reduced  frequency  of  0.2,  and  a 
pitch  amplitude  of  0.5°  calculated  with  TSP  and  Lin.  TSP.  For  the  TSP  solution  the  peaks  on  the  upper  side 
induced  by  the  shock  are  a  little  more  downstream  than  for  the  Lin.  TSP  solution.  This  effect  can  be 

produced  by  the  asymmetric  shock  motion  about  the  steady  shock  position  in  the  TSP  solution. 

Fig.  12  shows  the  unsteady  pressure  distribution  for  the  same  flow  parameters  calculated  with  vise.  TSP 
and  vise.  Lin.  TSP.  The  pressure  distributions  on  the  upper  and  lower  surface  are  well  comparable.  The 
peaks  of  the  pressure  distribution  in  the  vicinity  of  the  shock  on  the  upper  surface  calculated  with  vise. 
Lin.  TSP  are  a  little  higher  than  of  the  pressure  distribution  calculated  with  vise.  TSP.  This  effect  is 
valid  for  all  comparisons  for  pressure  distributions,  calculated  in  the  above  mentioned  way. 

Fig.  13  shows  the  unsteady  aerodynamic  coefficients  for  M=0.756  and  a =1.3°  and  a  range  of  reduced 
frequencies  between  0.05  and  0.5.  A  pure  elliptic  solution  of  the  TSP  code  for  the  flat  plate  is  indicated 
in  addition  to  the  above  mentioned  methods.  These  coefficients  are  also  derived  from  the  first  harmonic  of 
the  unsteady  pressure  distributions.  The  upper  part  of  the  figure  shows  the  coefficients  for  the  pitch 
motion  and  the  lower  part  those  for  the  heave  motion.  Because  of  the  neglect  of  the  -term  in  the  TSP 

equation  there  exists  a  small  discrepancy  for  "higher"  reduced  frequencies  between  D.L.  and  flat  plate 
results.  The  largest  difference  exists  beween  D.L.  and  TSP  results.  With  the  exception  of  the  imaginary 
part  of  Cma  and  the  real  part  of  Cmh  the  magnitude  of  the  difference  between  the  different  results  and  D.L. 
decreases  in  the  sequence  TSP,  Lin.  TSP,  vise.  Lin.  TSP  and  vise.  TSP.  The  most  striking  effect  is  the 
large  difference  between  TSP  and  vise.  TSP.  This  difference  is  on  the  average  larger  than  between  D.L.  and 
vise.  TSP.  This  effect  is  already  indicated  in  the  unsteady  pressure  distributions. 

Fig.  14,  15,  and  16  show  the  unsteady  aerodynamic  coefficients  against  frequency  for  Mach  numbers  0.76, 
0.78,  and  0.80  respectively.  Each  of  the  figures  shows  the  values  for  three  angles  of  attack,  calculated 
with  TSP;  plus  a  fourth  curve  which  represents  the  result  of  vise.  TSP  for  the  highest  incidence. 

Generally  the  values  for  the  real  and  imaginary  part  increase  with  increasing  incidence  except  for  a  =1.3° 
for  M=0. 78.  For  M=0.8  this  tendency  is  not  as  clear  as  for  M  =  0.765  and  M=0.78.  But  the  TSP  code  tends 
to  become  unstable  with  increasing  incidences.  The  aerodynamic  coefficients  calculated  with  vise.  TSP 
agrees  best  with  the  TSP  solution  for  the  lowest  incidence,  confirming  the  tendency  of  the  boundary  layer 
to  reduce  the  effective  incidence. 


Flutter  Calculation  Results 


To  show  the  influence  on  flutter  of  the  different  aerodynamic  input,  a  two-dimensional  dynamic  model  was 
selected.  Fig.  17  shows  the  mode  shapes  and  the  eigenfrequencies  of  this  model  which  are  chosen  in  such  a 
way,  that  they  represent  the  streamwise  motion  of  a  strip  in  the  outer  wing  section  of  a  backward  swept 
wing  of  a  transport  aircraft  for  the  two  eigenmodes  fundamental  bending  and  engine  pitch.  Mode  1  is 
characterized  by  a  large  heave  and  a  small  pitch  motion  with  a  nodal  point  upstream  of  the  section,  x*  is 
the  distance  between  the  nodal  point  and  the  1/4  point  of  the  section.  Mode  2  is  characterized  by  a  pitch 
motion  around  a  point  with  a  distance  x*  downstream  of  the  1/4  point  of  the  section. 
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The  flutter  results  were  obtained  for  a  variable  mass  ratio  corresponding  to  the  condition  in  a 
pressurized  windtunnel  for  selected  M-numbers  at  constant  stagnation  temperature.  Therefore  the  flutter 
results,  frequencies  and  damping  of  the  different  modes,  depend  on  the  stagnation  pressure  for  a  constant 
M-number.  The  most  important  coefficients  for  flutter  of  this  representative  section  are  Cla  and  Cl^  as 
long  as  the  nodal  point  of  the  second  mode  is  remote  from  the  1/4  point  of  the  section.  If  the  nodal  point 
of  the  second  mode  is  close  to  the  1/4  point  the  moment  coefficients  can  also  become  important  for  flutter. 

For  all  the  cases  the  mode  shapes  were  chosen  in  the  following  way:  Mode  1:  x*  =-10.;  Mode  2:  x*  =0.5. 

The  flutter  calculation  presented  in  Fig.  18  was  performed  with  a  chord  length  of  0.3  m  and  in  Fig.  19  with 

a  chord  length  of  0.6  m.  That  means  the  reduced  frequencies  used  in  Fig.  18  are  approximately  half  the 

reduced  frequencies  in  Fig.  19. 

For  the  first  case  (c=0.3  m),  the  critical  stagnation  pressures  that  were  calculated  for  various  Mach 
numbers  and  incidences  are  shown  on  Fig.  18  plotted  against  Mach  number.  The  Mach  number  ranges  from  0.756 
to  0.8,  the  incidences  from  0°  to  1.3°,  plus  a  case  for  M=0.82  and  a=-l°.  The  methods  that  were  used  to 
calculate  these  values  were:  D.L.,  TSP,  vise.  TSP,  lin.  TSP,  vise.  lin.  TSP. 

The  Mach  number  was  limited  at  the  lower  boundary  by  the  design  Mach  number,  and  at  the  upper  end  by  the 
breakdown  of  the  TSP  code.  The  upper  end  of  the  incidence  range  was  limited  by  either  the  breakdown  of  the 
TSP  code  or  the  onset  of  one-degree-of -freedom  flutter  at  zero  stagnation  pressure.  The 
one-degree-of-freedom  flutter  is  induced  for  this  small  reduced  frequency  by  the  magnitude  of  the  unsteady 
aerodynamic  coefficients:  the  imaginary  part  of  Cla  and  the  imaginary  part  of  Cmh  .  In  Fig.  18  two  flutter 
cases  are  involved:  one-degree-of-freedom  flutter  for  the  bending  mode  with  low  critical  stagnation 
pressures,  and  bending-torsion  flutter.  The  flutter  points  calculated  with  the  D.  L.  and  vise.  TSP  codes 
belong  to  the  bending-torsion  flutter  type,  but  the  flutter  points  calculated  with  the  TSP  code  are 
strongly  influenced  by  the  one-degree-of-freedom  flutter. 

Therefore  at  low  reduced  frequencies  the  noticeable  decrease  in  critical  stagnation  pressures, 
calculated  with  the  inviscid  codes,  is  caused  mainly  by  a  change  in  the  flutter  type.  Nevertheless  the 
curves  drawn  in  Fig.  18  for  constant  incidence  indicate  the  transonic  dip.  What  is  striking  here  is  the 
large  discrepancy  between  the  D.L.  solution  and  the  TSP  solutions  dependent  on  incidence,  compared  to  the 

small  differences  between  D.L.  and  vise.  TSP,  which  are  nearly  independent  of  incidence.  The  flutter  point 

calculated  with  TSP  for  M=0.8  and  a  =0.75  with  a  higher  stagnation  pressure  as  for  D.L.  is  not 
understandable,  but  not  a  failed  calculation. 

The  calculations  for  the  second  case  (c=0.6  m)  were  performed  with  the  same  aerodynamic  input.  For  this 
higher  reduced  frequency  range  no  one-degree-of-freedom  flutter  influence  appears,  as  in  Fig.  18. 

Therefore  the  limit  for  the  upper  M-number  and  incidence  boundary  was  defined  by  the  divergence  of  the  TSP 
code.  In  Fig.  19  a  typical  transonic  dip  is  visible  for  the  critical  pressure  points  calculated  with  the 
TSP  code,  which  is  also  due  to  the  dependence  on  incidence.  The  critical  pressure  points  calculated  with 
the  Lin.  TSP  with  and  without  the  steady  boundary  layer  also  show  the  dip  effect  in  the  expected  way.  But 

the  solutions  with  vise.  TSP  show  higher  values  for  the  critical  stagnation  pressure  for  the  M-numbers  0.78 

and  0.8  than  calculated  with  D.L.  This  result  can  be  explained  by  the  different  behaviour  of  the  pitch 
frequency  vs.  stagnation  pressure  for  vise.  TSP.  It  is  uncertain  whether  this  effect  is  realistic  since 
the  unsteady  terms  in  the  boundary  layer  were  neglected. 


Flutter  Simulation  Results 


Fig.  20  shows  the  typical  phases  of  the  flutter  simulation  approach.  In  the  first  part  of  the  steady 
phase,  incidence  and  profile  thickness  are  increased  from  zero  to  full  value,  to  start  the  calculation. 

Then  the  airfoil  is  held  motionless  until  a  steady  state  is  achieved.  In  the  second  phase  a  forced 
excitation  in  pitch  and  heave  is  applied: 

a  ( t )  =  a  sin  coi 
h(t)=  hsin(cvt+^) 

Frequency,  phase  and  ratio  of  amplitudes  of  the  excitation  are  the  same  as  those  of  the  flutter  mode  of 
the  flutter  calculation.  After  three  cycles  the  excitation  is  switched  off  leaving  the  section  in  free 
vibration.  In  order  to  keep  the  profile  in  a  mean  position  which  corresponds  to  the  steady  incidence,  an 
initial  stress  is  added  during  this  phase  to  compensate  the  static  unbalance.  Fig.  20  compares  a  TSP  to  a 
vise.  TSP  simulation  at  design  conditions.  The  time  responses  of  the  TSP  simulation  show  a  slightly 
unstable  development  while  the  vise.  TSP  responses  are  well  damped. 

In  Fig.  21  TSP  flutter  simulation  results  with  three  different  amplitudes  are  presented  for  M=0.756  and 
a  =1.3°.  The  left  side  of  the  figure  shows  the  time  responses  for  a  nearly  critical  stagnation  pressure, 
while  on  the  right  hand  side  a  supercritical  case  is  presented.  The  amplitude  growth  rates  become  larger 
for  higher  initial  excitation  amplitudes.  The  higher  harmonics  are  visible  only  in  the  time  histories  of 
the  moment  coefficient. 

Fig.  22  and  Fig.  23  show  the  frequency  and  damping  curves  vs.  stagnation  pressure  at  design  conditions 
for  the  flutter  model  with  0.3  m  and  0.6  m  respectively.  Flutter  calculation  results  are  compared  to 
flutter  simulation  results  for  three  different  amplitudes.  The  development  of  the  simulation  results  agree 
well  with  the  calculation  results.  The  simulation  approach  produces  flutter  at  about  5-10%  lower 
stagnation  pressures  depending  on  the  frequency  range  and  the  amplitude  ratio.  For  the  0.3  m  case  the 
flutter  pressure  reduces  with  increasing  amplitudes  while  for  the  0.6  m  case  the  dependence  of  amplitude  is 
low. 
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CONCLUSION 

For  moderate  transonic  Mach  numbers  and  incidences  the  steady  and  unsteady  pressure  distribution 
calculated  with  TSP  generally  agree  well  with  the  corresponding  ones  calculated  with  the  Euler  code.  The 
stability  of  the  TSP  code,  however,  is  weak. 

While  the  steady  flow  solutions  are  strongly  influenced  by  the  boundary  layer,  the  unsteady  solutions 
seem  to  be  affected  even  more  by  the  boundary  layer.  Inviscid  solutions,  which  are  matched  to  viscid 
steady  pressure  distributions  deviate  considerably  in  their  unsteady  behaviour.  Therefore  flutter 
calculations  that  use  airloads  calculated  by  nonlinear  methods  must  take  account  of  boundary  layer  effects, 
at  least  in  the  case  of  supercritical  profiles  or  profiles  with  strong  nonlinear  variation  of  Cl  with 
incidence,  or  if  the  boundary  layer  linearizes  the  Cl  vs.  incidence  variation  considerably. 

The  influence  of  higher  harmonics  in  the  unsteady  pressure  distribution,  and  the  nonlinear  dependence  on 
the  amplitude,  have  a  smaller  influence  on  flutter,  as  shown  by  the  comparison  between  the  results  of 
flutter  calculation  and  flutter  simulation. 


REFERENCES 

(11  Farmer,  M.G.,  Hanson,  P.W. 

Comparison  of  supercritical  and  conventional  wing  flutter  characteristics 
AIAA/ASME/SAE  17th  Structure, 

Structural  Dynamics  and  Material  Conference 
May  1976 

[21  Ballhaus,  KW.F.,  Bridgeman,  J.O. 

Numerical  Solution  Techniques  for  Unsteady  Transonic  Aerodynamics  Problems 
AGARD  Rep.  679,  March  1980 

[31  Tijdeman,  H. 

Investigations  of  the  transonic  flow  araound  oscillating  Airfoils 
NLR  TR  77090  U 

[41  Zimmermann,  H., 

Application  of  Transonic  Unsteady  Methods  for  Calculation  of  Flutter  Airloads 
AGARD,  Toulouse  1984 

[51  Couston,  M.,  Angelini,  J.-J.,  Mulak,  P. 

Application  de  1 'equation  des  petites  perturbations  transoniques  aux  calculs  d’ecoulements 

bidimensionel s  instationaires 

La  Recherche  Aerospatiale  No.  1979-5,  pp.  325-340 

[61  Thiede,  P.,  Elsholz,  E.,  Dargel ,  G. 

Rationelle  Berechnung  partieller  Ablosegebiete  auf  der  Basis  des  Grenzschichtkonzepts 

RiiFo  IV,  T/RF41  /80050/81449,  1981 

[7]  Zimmermann,  H.,  Vogel,  S. 

Application  of  Transonic  Unsteady  Methods  for  Calculation  of  Flutter  Airloads 
Second  International  Symposium  on  Aeroelasticity  and  Structural  Dynamics 
Aachen  1985 

[83  Dau,  K.,  Schulze,  B.,  Zimmermann,  H., 

Instationare  Luftkrafte  im  transsonischen  Bereich  und  ihre  Auswirkungen  auf  das  Flatterverhalten, 
Teil  III 
BMVg-FBWT  86-1 

[91  Vogel,  S. 

Ein  Beitrag  zur  Losung  der  FI attergleichung  unter  BerUcksichtigung  von  Servorsteuerung  und 
Flugregler 

ZFW  Band  1,  Heft  2,  1977 
[101  Bathe,  K.-J. 

Finite  Element  Procedures  in  Engineering  Analysis 
Prentice-Hall,  Inc.,  Englewood  Cliffs,  New  Jersey  07632,  1982 

[111  Jacob,  H.G. 

Rechnergestutzte  Optimierung  statischer  und  dynamischer  Systeme 
Springer-Verlag,  Berlin,  Heidelberg,  New  York,  1982 


2-7 


Relaxation  formulas:  6n+1  =  6°  ^'riv's 


Fig.  1  Direct  viscid-inviscid 
coupling  procedure 


Fig.  2  Indirect  viscid-inviscid 
coupling  procedure 


m 


Fig.  3  MBB-A3  profile 


ITER.  80-95 


Fig.  4  Iterative  development  of  the  boundary  layer  displacement  thickness 
MBB-A3,  M=0.8,  a=-0.2°  and  Re=17.7  106 
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Fig.  5  Steady  and  unsteady  pressure  distribution  for  MBB-A3  profile 
M=0.8,  <x=— 0.2°  and  k=0.104  (  pitch  motion  ) 


Fig.  6  Steady  pressure  distribution 
MBB-A3  (  design  conditions  ) 


Fig.  7  Steady  pressure  distribution 
MBB-A3  (  M=0.78  ) 
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Fig.  8  Unsteady  pressures,  (  pitch  motion  )  Fig.  9  Unsteady  pressures,  (  pitch  motion  ) 

calculated  with  TSP  calculated  with  vise.  TSP 

M=0.756,  a=1.3°  M=0.756,  a=1.3°  and  Re=  6.  106 


SURFACE  LOWER  SURFACE  UPPER  SURFACE  LOWER  SURFACE  UPPER  SURFACE  LOWER  SURFACE 
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Fig.  13  Unsteady  aerodynamic  coefficients  for 
M  =  0.756  and  a=1.3°  (Re=6.  106) 


JSP .  a  =  0.0°  •  TSP .  a  =  0.0°  •  TSP .  a  =  0.0° 

TSP .  a  =  0.5°  ©  TSP .  a  =  0.75°  ©  TSP .  a  =  0.5° 

TSP .  a  =  0.75°  *  TSP .  a  =  1.3°  *  TSP .  a  =  0.75 

VISC.  TSP.  a  =  0.75°  *  VISC.  TSP.  a  =  1.3°  *  VISC.  TSP.  a  =  0.75 
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Fig.  14  Unsteady  aerodynamic  coefficients  Fig.  15  Unsteady  aerodynamic  coefficients  Fig.  16  Unsteady  aerodynamic  coefficients 
M=0.765and  a=0° ;  0.5° ;  0.75°  M=0.78  and  a=0° ;  0.75° ;  1.3°  M=0.8  and  a=0° ;  0.5° ;  0.75° 


Pt  [bar] 
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MODE  1 


MODE  2 


Fig.  17  Dynamic  model  of  o  typical  section 
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FLUTTER  CALCULATION 
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Fig.  18  Flutter  results 
(chord  =  0.3m) 
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FLUTTER  CALCULATION 
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Fig.  19  Flutter  results  (chord  =  0.6m) 


Fig.  20  Flutter  simulation  results 

calculated  with  TSP  and  vise.  TSP 
M=0.756,a=1.3°,  Pt  =  0.55  bor 
(chord  =  0.3m) 


Fig.  21  Flutter  simulation  results  (chord  =  0.3m)  M=0.756,  a=1.3,° 
calculated  with  TSP 
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Fig.  22  Comparison  between  flutter  calculation 
and  flutter  simulation  results 
(chord  =  0.3m)  M=0.756,  a-  1.3° 
calculated  with  TSP 


•=  FLUTTER  CALCULATION  1 .MODE 
o=  FLUTTER  CALCULATION  2. MODE 
a=  FLUTTER  SIMULATION  a=0.1° 
x  =  FLUTTER  SIMULATION  a=0 . 5° 
o  =  FLUTTER  SIMULATION  5=1.0 


Fig.  23  Comparison  between  flutter  calculation 
and  flutter  simulation  results 
(chord  =  0.6m)  M=0.756,  a—  1.3° 
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o —  FLUTTER  CALCULATION  2. MODE 
A=  FLUTTER  SIMULATION  5=0.1° 
x=  FLUTTER  SIMULATION  5=0.5° 
□  =  FLUTTER  SIMULATION  5=1.0° 
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